/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Source_div
*
* %Identification
* Written by: Erik Knudsen 
* Date: November 11, 2009
* Origin: Risoe
* Release: McXtrace 0.1
*
* X-ray source with Gaussian or uniform divergence
*
* %Description
* A flat rectangular surface source with uniform or Gaussian divergence profile and focussing.
* If the parametere gauss is not set (the default) the divergence profile is flat
* in the range [-focus_ax,focus_ay]. If gauss is set, the focux_ax,focus_ay is considered
* the standard deviation of the gaussian profile.
* Currently focussing is only active for flat profile. The "focus window" is defined by focus_xw,focus_yh and dist.
* The spectral intensity profile is uniformly distributed in the energy interval defined by e0+-dE/2 or 
* by wavelength lambda0+-dlambda/2
* 
* Example: Source_div(xwidth=0.1, yheight=0.1, focus_aw=2, focus_ah=2, E0=14, dE=2, gauss=0)
*
* %Parameters
* xwidth:   [m]   Width of source.
* yheight:  [m]   Height of source.
* focus_aw: [rad] Standard deviation (Gaussian) or maximal (uniform) horz. width divergence.
* focus_ah: [rad] Standard deviation (Gaussian) or maximal (uniform) vert. height divergence.
* focus_xw: [m]   Width of sampling window
* focus_yh: [m]   Height of sampling window
* dist:     [m]   Downstream distance to place sampling target window
* E0:       [keV] Mean energy of X-rays.
* dE:       [keV] Energy half spread of X-rays. If gauss==0 dE is the half-spread, i.e. E\in[E0-dE,E0+dE], if gauss!=0 it's interpreted as the standard dev. 
* lambda0:  [AA]  Mean wavelength of X-rays (only relevant for E0=0).
* dlambda:  [AA]  Wavelength half spread of X-rays.
* gauss:    [1]   Criterion: 0: uniform, 1: Gaussian distribution of energy/wavelength.
* gauss_a:  [1]   Criterion: 0: uniform, 1: Gaussian divergence distribution.
* flux:     [1/(s * mm**2 *mrad**2 * energy unit)] flux per energy unit, Angs or keV.
* randomphase: [1] If !=0 the photon phase is chosen randomly.
* phase:    [1]   Value of the photon phase (if randomphase==0).
* spectrum_file: [string] File from which to read the spectral intensity profile
* focus_ar: [rad] Standard deviation (Gaussian) or maximal (uniform) radial divergence.
* radius:   [m]   Radius of circular source
* verbose:  [0/1] Show more information
*
* %End
*******************************************************************************/

DEFINE COMPONENT Source_div

SETTING PARAMETERS (string spectrum_file="NULL", xwidth=0, yheight=0, dist=0,
    focus_xw=0, focus_yh=0,focus_aw=0, focus_ah=0, focus_ar=0, radius=0,
    E0=0, dE=0, lambda0=0, dlambda=0, flux=0, gauss=0, gauss_a=0, int randomphase=1, phase=0, int verbose=0)

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 

SHARE
%{
  %include "read_table-lib"
%}

DECLARE
%{
  double p_init;
  double K;
  double dK;
  double xmin;
  double xmax;
  double xw_2;
  double focus_xw_2;
  double ymin;
  double ymax;
  double yh_2;
  double focus_yh_2;
  double pmul;
  double pint;
  t_Table T;
  int spectrum_from_file;
%}

INITIALIZE
%{
  focus_xw_2=focus_xw/2.0;
  focus_yh_2=focus_yh/2.0;
  xmin=-xwidth/2.0;
  ymin=-yheight/2.0;
  xmax=xwidth/2.0;
  ymax=yheight/2.0;

  if(radius==0 && ( xmax==xmin && ymin==xmax) ){
    fprintf(stderr,"ERROR (%s): No meaningful source area set. Either set radius or xwidth and yheight.\n",NAME_CURRENT_COMP);
    exit(-1);
  }

  /*flag if we are using a datafile*/
  spectrum_from_file=(spectrum_file && strcmp(spectrum_file,"NULL") && strlen(spectrum_file));

  if (spectrum_from_file){
    /*read spectrum from file*/
    int status=0;
    if ( (status=Table_Read(&(T),spectrum_file,0))==-1){
      fprintf(stderr,"ERROR (%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,spectrum_file?spectrum_file:"");
      exit(-1);
    }
    /*data is now in table t*/
    /*integrate to get total flux, assuming raw numbers have been corrected for measuring aperture*/
    int i;
    pint=0;
    for (i=0;i<T.rows-1;i++){
      pint+=((T.data[i*T.columns+1]+T.data[(i+1)*T.columns+1])/2.0)*(T.data[(i+1)*T.columns]-T.data[i*T.columns]);
    }
    if (verbose){
      printf("INFO (%s): Integrated intensity radiated is %g pht/s\n",NAME_CURRENT_COMP,pint);
      if(E0) printf("INFO (%s):, E0!=0 -> assuming intensity spectrum is parametrized by energy [keV]\n",NAME_CURRENT_COMP);
    }
  }else if (!E0 && !lambda0){
    fprintf(stderr,"ERROR (%s): Error: Must specify either wavelength or energy distribution\n",NAME_CURRENT_COMP);
    exit(-1);
  }  

  /*calculate the X-ray weight from the flux*/
  if (flux){
    pmul=flux;
  }else{
    pmul=1;
  }
  pmul*=1.0/((double) mcget_ncount());

  if( dist==0 && ( focus_xw!=0 || focus_yh!=0 )){
    fprintf(stderr,"ERROR (%s): Cannot have focus sampling window (focus_xw x focus_yh) = (%g x %g) with dist=0.\n",NAME_CURRENT_COMP);
    exit(-1);
  }

  /*check if divergence limits are compatible with focus_xw, focus_yh*/
  double maxdivh,maxdivv;
  if(focus_xw!=0){
    maxdivh=atan((xwidth+focus_xw)/dist);
    if (focus_aw>maxdivh){
      focus_aw=maxdivh;
      if (verbose){
        fprintf(stderr,"WARNING (%s): sampling width does not support full divergence. Adjusting to focus_aw=%g rad\n",NAME_CURRENT_COMP,focus_aw);
      }
    }
  }
  if(focus_yh!=0){
    printf("I am here\n");
    maxdivv=atan((yheight+focus_yh)/dist);
    if (focus_ah>maxdivv){
      focus_ah=maxdivv;
      if (verbose){
        fprintf(stderr,"WARNING (%s): sampling height does not support full divergence. Adjusting to focus_ah=%g rad\n",NAME_CURRENT_COMP,focus_ah);
      }
    }
  }
%}

TRACE
%{
  double kk,theta_x,theta_y,l,e,k;
  p=pmul;

  if (spectrum_from_file){
    double pp=0;
    while (pp<=0){ 
      l=T.data[0]+ (T.data[(T.rows-1)*T.columns] -T.data[0])*rand01();
      pp=Table_Value(T,l,1);
    }
    p*=pp;
    /*if E0!=0 the tabled value is assumed to be energy in keV*/
    if (E0!=0){
      k=E2K*l;
    }else{
      k=(2*M_PI/l);
    }
  }else if (E0){
    if(!dE){
      e=E0;
    }else if (gauss){
      e=E0+dE*randnorm();
    }else{
      e=randpm1()*dE + E0;
    }
    k=E2K*e;
  }else if (lambda0){
    if (!dlambda){
      l=lambda0;
    }else if (gauss){
      l=lambda0+dlambda*randnorm();
    }else{
      l=randpm1()*dlambda*0.5 + lambda0;
    }
    k=(2*M_PI/l);
  }
  /*pick a point of origin*/
  if(!radius){
    x=xmin+rand01()*xwidth;
    y=ymin+rand01()*yheight;
    z=0;
  }else{
    double r=radius*sqrt(rand01());
    double th=rand01()*2*M_PI;
    x=r*cos(th);
    y=r*sin(th);
    z=0;
  }
  /*pick a direction*/
  if (focus_aw!=0 || focus_ah!=0){
    if (!gauss_a){
      /*find limits of uniform sampling scheme for vertical divergence.
        thetav should be acos(1-2*U) for U\in[0,1]. for theta measured from vertical axis
        we only use a sub-interval for U and measure from horizontal plane.*/
      double sample_lim1,u2;
      sample_lim1=(1-cos(M_PI_2 - focus_ah/2.0))*0.5;
      u2=randpm1()*(sample_lim1-0.5) + 0.5;
      theta_x = randpm1()*focus_aw/2.0;
      theta_y = acos(1-2*u2) - M_PI_2;
    }else{
      theta_x=randnorm()*focus_aw;
      theta_y=randnorm()*focus_ah;
    }

    kx=tan(theta_x);
    ky=tan(theta_y);
    kz=1.0;
    NORM(kx,ky,kz);
    kz*=k;
    kx*=k;
    ky*=k;
  }else if (focus_ar>0){
    /*radial divergence profile*/
    double theta,psi,kp;
    if(!gauss_a){
      double sample_lim1,u2;
      sample_lim1=(1-cos(focus_ar/2.0))*0.5;
      u2=rand01()*(sample_lim1);
      psi=acos(1-2*u2);
    }else{
      do{
        psi=abs(randnorm()*focus_ar);
      } while (psi>M_PI);
    }
    theta=rand01()*2*M_PI;

    kz=k*cos(psi);
    kp=sqrt(k*k-kz*kz);
    kx=kp*sin(theta);
    ky=kp*cos(theta);
  }else{
    /*0 dviergence - source is mathematically collimated*/
    kz=k;
    kx=ky=0;
  }

  /*set polarization and phase.*/
  Ex=0;Ey=0;Ez=0;
  if (!randomphase){
    phi=phase;
  }else{
    phi=rand01()*M_2_PI;
  }
%}

MCDISPLAY
%{
  double dist_display=1;
  if (dist_display>dist){
    dist_display=dist;
  }
  multiline(5, -xwidth/2.0, -yheight/2.0, 0.0,
                xwidth/2.0, -yheight/2.0, 0.0,
                xwidth/2.0,  yheight/2.0, 0.0,
               -xwidth/2.0,  yheight/2.0, 0.0,
               -xwidth/2.0, -yheight/2.0, 0.0);
  if (focus_aw){
    dashed_line(0,0,0, tan(focus_aw/2.0)*dist_display,0,dist_display,4);
    dashed_line(0,0,0,-tan(focus_aw/2.0)*dist_display,0,dist_display,4);
  }
  if (focus_ah){
    dashed_line(0,0,0,0, tan(focus_ah/2.0)*dist_display,dist_display,4);
    dashed_line(0,0,0,0,-tan(focus_ah/2.0)*dist_display,dist_display,4);
  }
%}

END