/*
 * passage.c
 *
 * Written by Ricardo Yanez <ricardo.yanez@calel.org>
 *
 */

#include <stdio.h>
#include <math.h>
#include <range.h>
#include <stdlib.h>

/* Calculate energy loss of alpha particles in scint */

int main (int argc, char** argv) {

  int i;
  double dens;
  double ptot, mp;
  double theta;
  double emin, emax, ecurrent, eincr;
  double ein;
  double eout, pout, p_target_0;
  double thickness, thickness_0, error;

  /* density of air at 1 atm (see README file) */
  dens = 1.03; /* g/cm^3 */
  dens = dens*1000.0; /*  mg/cm^3 */

  /* define absorber (scint) */
  nelem = 2;
  absorb[0].z = 1;  absorb[0].a =  1; absorb[0].w = 2*1; /* H */
  absorb[1].z = 6;  absorb[1].a = 12; absorb[1].w = 1*12; /* C */

  /* If you want to be more accurate you can use molecular weights instead
     of mass numbers when defining the weights */

  mp = 938.27;
  thickness_0 = atof(argv[1]); // thickness in gm/cm^2
  p_target_0 = atof(argv[2]); // minimum trackable momentum in MeV/c
  for (theta = 1.0; theta < 90.; theta += 1.0) {
    emin = 0.0;
    emax = 1000.;
    ecurrent = (emax + emin)/2.0;
    eincr = (emax - emin)/2.0;
    for (i = 0 ; i < 20 ; i++) {
      ein = ecurrent;
      ptot = sqrt((mp + ein)*(mp + ein) - mp*mp);
      thickness = thickness_0/sin(theta*3.14159/180.);
      eout = passage(0,1,1,-1,0,0,ein, thickness, &error);
      pout = sqrt((mp + eout)*(mp + eout) - mp*mp);
      //printf("ein = %6.11f, ptot = %6.11f, eout = %6.11f, pout = %6.11f\n", ein, ptot, eout, pout);
      eincr = eincr/2.0;
      if (pout > p_target_0) {ecurrent -= eincr;} else {ecurrent += eincr;}
    }

    //printf("theta = %6.2f, ptot = %6.2f\n", theta, ptot);
    printf("%6.2f %6.2f\n", theta, ptot);

  }

  return 0;

}