#include <iostream>
#include <fstream>
#include <math.h>
#include <cstdlib>

using namespace std;


double c = 3e8; // [m/s] predkosc swiatla
double R = 0.8875e-15; // [m] promien protonu
double m = 1.67262178e-27; // [kg] masa protonu
double e = 1.602176565e-19; // [C] ladunek elementarny
double k = 8.987e9; // [N m^2 C^-2] stala Coulomba


// Zaleznosc kata rozpraszania od parametru zderzenia dla kulek bilardowych
double rozpraszanie1(double b){
  double fi;
  if( fabs(b) > 2*R ) 
    fi = 0.;
  else 
    fi = 180. - 360.*asin(b/(2*R))/3.1415;
  return fi;
}


// Zaleznosc kata rozpraszania od parametru zderzenia dla rozpraszania sila
// Coulomba
double rozpraszanie2(double b){
  double A = k*e*e/m;
  double x0 = b;
  double y0 = -1000*R;
  double vx0 = 0;
  double vy0 = c;

  double r, x, y, vx, vy, ax, ay, fi, dt;
  x = b;
  y = y0;
  vx = vx0;
  vy = vy0;

  // Powtarzamy, az pocisk oddali sie na taka odleglosc, jak byla na poczatku
  while (x*x + y*y <= x0*x0 + y0*y0) {
    r = sqrt(x*x + y*y); // odleglosc od tarczy
    ax = A*x/(r*r*r); // przyspieszenie w kierunku x
    ay = A*y/(r*r*r); // przyspieszenie w kierunku y
    dt = 100*r*r; // przedzial czasu - tym mniejszy im wieksza sila
    vx = ax*dt + vx; // zmiana predkosci
    vy = ay*dt + vy; // zmiana predkosci
    x = 0.5*ax*dt*dt + vx*dt + x; // zmiana polozenia
    y = 0.5*ay*dt*dt + vy*dt + y; // zmiana polozenia
  }

  // wyliczenie kata rozproszenia
  fi = fabs(atan2(vx, vy));
  return fi*180./3.14159265;
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// Czesc pierwsza cwiczenia - zaleznosc kata rozpraszania od parametru zderzenia
void czesc1(){
  ofstream plik("zaleznosc.xls");
  for(double x = 0; x<5; x = x + 0.01){
    double b = R*x;
    double fi1 = rozpraszanie1(b);
    double fi2 = rozpraszanie2(b);
    plik << x << "\t" << fi1 << "\t" << fi2 << endl;
    cout << x << "\t" << fi1 << "\t" << fi2 << endl;
  }
  plik.close();
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////



// Losowanie liczby z podanego przedzialu
double losuj(double xmin, double xmax){
  double zakres = xmax - xmin;
  double x = xmin + zakres * rand()/RAND_MAX;
  return x;
}


// Definicje histogramu
const int N = 20; // liczba binow
const double vmin = 0.; // minimalna wartosc
const double vmax = 180; // maksymalna wartosc
const double zakres = vmax - vmin; // dlugosc zakresu
int histogram[N] = {0}; // definicja histogramu i wypelnienie go zerami

// Funkcja wypelniajaca histogram
void wypelnij(double v){
  for (int i=0; i<N; i++) {
    if ( (v >= vmin + i*zakres/N) && (v < vmin + (i+1)*zakres/N) )
      histogram[i]++;
  }
}


// Wyswietla na ekran (konsole) zawartosc histogramu
void wyswietl(){
  for (int i=0; i<N; i++){
    cout << vmin + (i+0.5)*zakres/N << "\t" << histogram[i] << endl;
  }
}


// Zapisuje zawartosc histogramu do pliku o podanej nazwie
void zapisz(const char nazwa[]){
  ofstream plik(nazwa);
  for (int i=0; i<N; i++){
    plik << vmin + (i+0.5)*zakres/N << "\t" << histogram[i] << endl;
  }
  plik.close();
}


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// Czesc pierwsza cwiczenia - porownanie do wynikow LHC
void czesc2(){

  // Probujemy bardzo wiele losowych przypadkow rozpraszania
  for (double a=0; a<100000; a++)
  {    
    // Losujemy polozenie pocisku
    double xymax = 3*R;
    double x = losuj(-xymax, xymax);
    double y = losuj(-xymax, xymax);
    double b = sqrt(x*x + y*y); // Liczymy parametr zderzenia
    double fi = rozpraszanie2(b); // Liczymy kat rozpraszania
    wypelnij(fi); // Wypelniamy histogram
  }

  wyswietl(); // Wyswietlamy histogram na ekran
  zapisz("histogram.xls"); // Zapisujemy histogram do pliku
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////


// Program glowny
int main(){
  czesc1();
}