/* Aldiffr.c - Diffraction from Al foil on boundary elements */

/* This file is an addon to GPT to simulate electron diffraction in thin foils and was written by:
 * Pietro Musumeci - UCLA Physics Department
 *
 * Example use:
 * aldiffr("wcs","I","scat",0.0);
 * scatterplate("wcs","z",1, 0.002,0.002) scatter="scat";
 *
 */

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include "elem.h"


struct aldiffr_info
{
  double P ;
  double nmin ;

  struct scatter_info scatinfo ;
} ;

void aldiffr_scat(gptpar *par,double t,double dt, gpttrajectory *trajectory,void *info) ;


void aldiffr_init(gptinit *init)
{
  struct aldiffr_info *info ;
  char *name ;
  int numarg ;

  gptbuildECS( init ) ;

  numarg = gptgetargnum(init) ;

  if( numarg!=2 && numarg!=3 )
    gpterror( "Syntax: %s(ECS,name,P,[nmin])\n", gptgetname(init) ) ;

  info = gptmalloc( sizeof(struct aldiffr_info) ) ;

  name       = gptgetargstring(init,1) ;
  info->P    = gptgetargdouble(init,2) ;

  if( numarg==3 )
    info->nmin = gptgetargdouble(init,3) ;
  else
    info->nmin = 0.0 ;

  if( info->P<0 || info->P>1 )
    gpterror( "Probability P must be between 0 and 1." ) ;
  if( info->nmin < 0 )
    gpterror( "Nmin must be nonnegative." ) ;

  /* Install all functions */
  gptscatterinit(init,&info->scatinfo,name) ;
  gptinstallscatterfnc(name,aldiffr_scat,info) ;
}


void aldiffr_scat(gptpar *par,double t,double dt, gpttrajectory *ptraj,void *vinfo)
{
  struct aldiffr_info *info = (struct aldiffr_info *)vinfo ;

  double Gnew, lenGBnew,angle,rn1, rn2, rn3, rn4, d, theta;
  double newdr[3], GBnew[3], vnew[3], rnew[3];
  int i,nscatt, j, sign;
  double scatterprob[6];
  double ringprob[6];
  double dring[5];

  /* how many scattering events? */
  /* Depends on the length of the foil and the cross section */
  /* The following is good for a 4 MeV beam into a 200 nm foil */

  scatterprob[0]=0.411;
  scatterprob[1]=0.777;
  scatterprob[2]=0.939;
  scatterprob[3]=0.987;
  scatterprob[4]=0.998;
  scatterprob[5]=1;

  /* Probability of scattering with a given angle */ 

  ringprob[0]=0;
  ringprob[1]=0.5;
  ringprob[2]=0.721;
  ringprob[3]=0.838;
  ringprob[4]=0.965;
  ringprob[5]=1;

  /* interplanar distance in Angstrom */

  dring[0]= 2.34;
  dring[1]= 2.02;
  dring[2]= 1.43;
  dring[3]= 1.22;
  dring[4]= 1.17;
 

  /* Generate new particles */

  Gnew = ptraj->Gint ;
  lenGBnew = gptVECLEN(ptraj->GBint) ;
  newdr[0] = ptraj->ndr[0];
  newdr[1] = ptraj->ndr[1]; 

  /* Determine number of scattering events */

  rn1 = dblpprand();

  if (rn1 < scatterprob[0]) return;
  for(i=1; i<6; i++) if(scatterprob[i-1] <= rn1 && rn1 < scatterprob[i]) nscatt = i;
  
  /* extract random phi angle */

  rn3 = dblpprand();
  angle = rn3*2*3.1415;

  for( j = 0 ; j < nscatt ; j++)  
  {
  
  /* Determine angle of scattering (theta) */
  rn2 = dblpprand();


  for(i=0; i<5; i++) if (ringprob[i]<= rn2 && rn2 < ringprob[i+1]) d = dring[i]; 

  theta = 0.024263/lenGBnew / d ;
  
  rn4 = dblpprand()-0.5;

  if (rn4 <= 0 ) sign = -1;
  if (rn4 > 0) sign = 1; 
 
  /* Forward scatter */
  newdr[0] = newdr[0]+theta*sign*cos(angle);	
  newdr[1] = newdr[1]+theta*sign*sin(angle);

  }

  newdr[2] = ptraj->ndr[2];
  for(i=0 ; i<3 ; i++) GBnew[i] = newdr[i] * lenGBnew ;  
  for(i=0 ; i<3 ; i++) vnew[i] = gpt_c*GBnew[i]/Gnew ;
  for(i=0 ; i<3 ; i++) rnew[i] = ptraj->P[i]+vnew[i]*dt*(1-ptraj->lambda) ;

  gptaddparmqn(gptgetparset("beam"), rnew, GBnew, par->m, par->q, par->n) ;

  /* Remove original particle and exit */
  gptremoveparticle(par) ;
  return;
}