/*
 * 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, range, lrange;
  double lrange_target_0, lrange_target;

  /* 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;
  lrange_target_0 = atof(argv[1]); // thickness in mm
  for (theta = 1.0; theta < 90.; theta += 1.0) {
    lrange_target = lrange_target_0/sin(theta/180.*3.14159);
    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);
      range = rangen(0,1,1,-1,0,0,ein);
      lrange = range/dens*10.0; /* mm */
      // printf("ein = %6.11f, ptot = %6.11f, range = %6.11f, lrange = %6.11f\n", ein, ptot, range, lrange);
      eincr = eincr/2.0;
      if (lrange > lrange_target) {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;

}