/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2006, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Source_gen
*
* %I
* Written by: Emmanuel Farhi, Kim Lefmann
* Date: Aug 27, 2001
* Version: $Revision: 1.26 $
* Origin: ILL/Risoe
* Release: McStas 1.9.1
* Modified by: EF, Aug 27, 2001 ; can use Energy/wavelength and I1
* Modified by: EF, Sep 18, 2001 ; corrected illumination bug. Add options
* Modified by: EF, Oct 30, 2001 ; minor changes. mccompcurname replaced
*
* Circular/squared neutron source with flat or Maxwellian energy/wavelength
* spectrum (possibly spatially gaussian)
*
* %D
* This routine is a neutron source (rectangular or circular), which aims at
* a square target centered at the beam (in order to improve MC-acceptance
* rate). The angular divergence is then given by the dimensions of the
* target.
* The neutron energy/wavelength is distributed between Emin=E0-dE and
* Emax=E0+dE or Lmin=Lambda0-dLambda and Lmax=Lambda0+dLambda. The I1 may
* be either arbitrary (I1=0), or specified in neutrons per steradian per
* square cm per Angstrom. A Maxwellian spectra may be selected if you give
* the source temperatures (up to 3). Finally, a file with the flux as a
* function of the wavelength [lambda(AA) flux(n/s/cm^2/st/AA)] may be used
* with the 'flux_file' parameter. Format is 2 columns free text.
* The source shape is defined by its radius, or can alternatively be squared
* if you specify non-zero h and w parameters.
* The beam is spatially uniform, but becomes gaussian if one of the source
* dimensions (radius, h or w) is negative, or you set the gaussian flag.
* Divergence profiles are triangular.
* The source may have a thickness, which will broaden the default zero time
* distribution.
* For the Maxwellian spectra, the generated intensity is dPhi/dLambda (n/s/AA)
* For flat spectra, the generated intensity is Phi (n/s).
*
* Usage example:
*   Source_gen(radius=0.1,Lambda0=2.36,dLambda=0.16, T1=20, T2=38, I1=1e13)
*   Source_gen(h=0.1,w=0.1,Emin=1,Emax=3, I1=1e13, verbose=1, gaussian=1)
*   EXTEND
*   %{
*      t = rand0max(1e-3); // set time from 0 to 1 ms for TOF instruments.
*   %}
*
* <b>Some sources parameters:</b>
* PSI cold source T1=150.42,I1=3.67e11 T2=38.74,I2=3.64e11
*                 T3=14.84,I3=0.95e11
* ILL VCS cold source T1=216.8,I1=1.24e+13,T2=33.9,I2=1.02e+13
*        (H1)         T3=16.7 ,I3=3.0423e+12
* ILL HCS cold source T1=413.5,I1=10.22e12,T2=145.8,I2=3.44e13
*        (H5)         T3=40.1 ,I3=2.78e13
* ILL Thermal tube    T1=683.7,I1=0.5874e+13,T2=257.7,I2=2.5099e+13
*        (H12)        T3=16.7 ,I3=1.0343e+12
* ILL Hot source      T1=1695,I1=1.74e13,T2=708,I2=3.9e12
*
* %VALIDATION
* Feb 2005: output cross-checked for 3 Maxwellians against VITESS source
*           I(lambda), I(hor_div), I(vert_div) identical in shape and absolute values
* Validated by: K. Lieutenant
*
* %P
* radius: [m]   Radius of circle in (x,y,0) plane where neutrons
*               are generated. You may also use 'h' and 'w' for a square source
* dist:   [m]   Distance to target along z axis.
* xw:     [m]   Width(x) of target.
*                 If dist=0.0, xw = horz. div  in deg
* yh:     [m]   Height(y) of target.
*                 If dist=0.0, yh = vert. div in deg
*
* Energy or wavelength range must be defined by giving min and max value or average and spread:
* Emin:     [meV] Minimum energy of neutrons
* Emax:     [meV] Maximum energy of neutrons
* E0:       [meV] Mean energy of neutrons.
* dE:       [meV] Energy spread of neutrons, half width.
* Lmin:     [AA]  Minimum wavelength of neutrons
* Lmax:     [AA]  Maximum wavelength of neutrons
* Lambda0:  [AA]  Mean wavelength of neutrons.
* dLambda:  [AA]  Wavelength spread of neutrons,half width
*
* Optional parameters:
* h:        [m]   Source y-height, then does not use radius parameter
* w:        [m]   Source x-width, then does not use radius parameter
* length:   [m]   Source z-length, not anymore flat
* T1:       [K]   Temperature of the Maxwellian source, 0=none
* I1:       [1/(cm**2*st*AA)] Source flux per solid angle, area and Angstrom
* T2:       [K]   Second Maxwellian source Temperature, 0=none
* I2:       [1/(cm**2*st*AA)] Second Maxwellian Source flux
* T3:       [K]   Third Maxwellian source Temperature, 0=none
* I3:       [1/(cm**2*st*AA)] Third Maxwellian Source flux
* flux_file:[str] name of a file that contains the wavelength distribution of
*                 the flux in [1/(s*cm**2*st*AA)]. Format is [lambda flux]
*                 if a file is given, temperature and intensity values are ignored
* gaussian: [0/1] Use gaussian divergence beam, normal distributions
* verbose:  [0/1] display info about the source. -1 unactivate source.
*
* Output parameters:
* jmax:     [1]   number of data points in flux distribution array
*
* %L
* <a href="http://www.ill.fr/Yellowbook">The ILL Yellow book</a>
* %L
* P. Ageron, Nucl. Inst. Meth. A 284 (1989) 197
* %E
******************************************************************************/


DEFINE COMPONENT Source_gen
DEFINITION PARAMETERS (flux_file=0)
SETTING PARAMETERS (radius=0.0, dist=57.296, xw=0, yh=0, E0=0, dE=0, Lambda0=0, dLambda=0, I1=0,
                    h=0, w=0, gaussian=0, verbose=0, T1=0,
                    Lmin=0,Lmax=0,Emin=0,Emax=0,T2=0,I2=0,T3=0,I3=0,length=0)
OUTPUT PARAMETERS (p_in,lambda0, lambda02, L2P,lambda_0,d_lambda,lambda0b, lambda02b, L2Pb,lambda0c, lambda02c, L2Pc, pTable)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

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

DECLARE
%{
 double p_in;
 double lambda_0=0,d_lambda=0;
 double lambda0, lambda02, L2P;     /* first Maxwellian source */
 double lambda0b, lambda02b, L2Pb;  /* second Maxwellian source */
 double lambda0c, lambda02c, L2Pc;  /* third Maxwellian source */
 t_Table pTable;
%}

INITIALIZE
%{
  double delta_lambda, source_area, k;

  if (radius == 0 && h == 0 && w == 0)
  {
    fprintf(stderr,"Source_gen: %s: Error: Please precise source geometry (radius, h, w)\n",
      NAME_CURRENT_COMP);
    exit(-1);
  }
  if (xw*yh == 0)
  {
    fprintf(stderr,"Source_gen: %s: Error: Please precise source target (xw, yh)\n",
      NAME_CURRENT_COMP);
    exit(-1);
  }
  if (Emin < 0 || Emax < 0 || Lmin < 0 || Lmax < 0 || E0 < 0 || dE < 0 || Lambda0 < 0 || dLambda < 0)
  {
    fprintf(stderr,"Source_gen: %s: Error: Negative average\n"
                   "            or range values for wavelength or energy encountered\n",
                   NAME_CURRENT_COMP);
    exit(-1);
  }
  if ((Emin == 0 && Emax > 0) || (dE > 0 && dE >= E0))
  {
    fprintf(stderr,"Source_gen: %s: Error: minimal energy cannot be less or equal zero\n",
      NAME_CURRENT_COMP);
    exit(-1);
  }

  if ((Emax >= Emin) && (Emin > 0))
  { E0 = (Emax+Emin)/2;
    dE = (Emax-Emin)/2;
  }
  if ((E0 > dE) && (dE >= 0))
  {
    Lmin = sqrt(81.81/(E0+dE)); /* Angstroem */
    Lmax = sqrt(81.81/(E0-dE));
  }
  if (Lmax > 0)
  { Lambda0 = (Lmax+Lmin)/2;
    dLambda = (Lmax-Lmin)/2;
  }
  if ((Lambda0 >= dLambda) && (dLambda > 0))
  { lambda_0 = Lambda0;
    d_lambda = dLambda;
  }
  else
  { fprintf(stderr,"Source_gen: %s: Error: Wavelength range %.3f +/- %.3f AA calculated \n",
      NAME_CURRENT_COMP, Lambda0, dLambda);
    fprintf(stderr,"- whole wavelength range must be >= 0 \n");
    fprintf(stderr,"- range must be > 0; otherwise intensity gets zero, use other sources in this case \n\n");
    exit(-1);
  }

  radius = fabs(radius); w=fabs(w); h=fabs(h);  I1=fabs(I1);
  lambda_0=fabs(lambda_0); d_lambda=fabs(d_lambda);
  xw = fabs(xw); yh=fabs(yh); dist=fabs(dist);

  if (dist == 0)
  {
    fprintf(stderr,"Source_gen: %s: warning: focusing distance is null.\n"
                   "            xw and yh interpreted as divergence in deg\n",
                   NAME_CURRENT_COMP);
  }

  delta_lambda = 2*d_lambda;
  Lmin = lambda_0 - d_lambda; /* Angstroem */
  Lmax = lambda_0 + d_lambda;

  if ((I1 > 0  && T1 > 0) || (flux_file && strlen(flux_file) > 0))
  {
    if ((h == 0) || (w == 0))
      source_area = radius*radius*PI*1e4; /* circular cm^2 */
    else
      source_area = h*w*1e4; /* square cm^2 */
    p_in = source_area*delta_lambda;
  }
  else
    p_in = 1.0/4/PI; /* Small angle approx. */
  p_in /= mcget_ncount();

  if (flux_file && strlen(flux_file) > 0) {
    if (Table_Read(&pTable, flux_file, 1) <= 0) /* read 1st block data from file into pTable */
      exit(fprintf(stderr, "Source_gen: %s: can not read file %s\n", NAME_CURRENT_COMP, flux_file));
    else { /* put table in Log scale */
      int i;
      if (pTable.columns != 2) exit(fprintf(stderr, "Source_gen: %s: Flux file %s should contain 2 columns\n", NAME_CURRENT_COMP, flux_file));
      for (i=0; i<pTable.rows; i++) {
        double val = Table_Index(pTable, i,1);
        if (val <= 0) fprintf(stderr, "Source_gen: %s: File %s has negative flux at row %i\n", NAME_CURRENT_COMP, flux_file, i+1);
        val = log(val > 0 ? val : 1/FLT_MAX);
        Table_SetElement(&pTable, i, 1, val);
      }
    }
  } else
  {
    k  = 1.38066e-23; /* k_B */
    if (T1 > 0)
    {
      lambda0  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T1);
      lambda02 = lambda0*lambda0;
      L2P      = 2*lambda02*lambda02;
    }
    else
      { lambda0 = lambda_0; }

    if (T2 > 0)
    {
      lambda0b  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T2);
      lambda02b = lambda0b*lambda0b;
      L2Pb      = 2*lambda02b*lambda02b;
    }
    else
      { lambda0b = lambda_0; }

    if (T3 > 0)
    {
      lambda0c  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T3);
      lambda02c = lambda0c*lambda0c;
      L2Pc      = 2*lambda02c*lambda02c;
    }
    else
      { lambda0c = lambda_0; }
  }

  if (verbose)
  {
    printf("Source_gen: component %s ", NAME_CURRENT_COMP);
    if ((h == 0) || (w == 0))
      printf("(disk)");
    else
      printf("(square)");
    printf("\n            spectra ");
    printf("%.3f to %.3f AA (%.3f to %.3f meV)", Lmin, Lmax, 81.81/Lmax/Lmax, 81.81/Lmin/Lmin);
    if (gaussian)
      printf(", gaussian divergence beam");
    printf("\n");
    if (flux_file && strlen(flux_file) > 0)
    { printf("  File %s for flux distribution used. Flux is dPhi/dLambda in [n/s/AA]. \n", flux_file);
      Table_Info(pTable);
    }
    else if (T1 && I1)
    { if (T1 != 0)
        printf("            T1=%.1f K (%.3f AA)", T1, lambda0);
      if (T2*I2 != 0)
        printf(", T2=%.1f K (%.3f AA)", T2, lambda0b);
      if (T3*I3 != 0)
        printf(", T3=%.1f K (%.3f AA)", T3, lambda0c);
      if (T1) printf("\n");
      printf("  Flux is dPhi/dLambda in [n/s/AA].\n");
    }
    else
    { printf("  Flux is Phi in [n/s].\n");
    }
  }
  else
    if (verbose == -1)
      printf("Source_gen: component %s unactivated", NAME_CURRENT_COMP);
%}

TRACE
%{
  double theta0,phi0,theta1,phi1,chi,theta,phi,v,r, lambda;
  double tan_h, tan_v, Maxwell, lambda2, lambda5;

  if (verbose >= 0)
  {

    z=0;

    if ((h == 0) || (w == 0))
    {
      chi=2*PI*rand01();                          /* Choose point on source */
      r=sqrt(rand01())*radius;                    /* with uniform distribution. */
      x=r*cos(chi);
      y=r*sin(chi);
    }
    else
    {
      x = w*randpm1()/2;   /* select point on source (uniform) */
      y = h*randpm1()/2;
    }
    if (length != 0)
      z = length*randpm1()/2;

    if (dist == 0)
    {
      theta0 =  DEG2RAD*xw;
      phi0   =  DEG2RAD*yh;
      theta1 = -DEG2RAD*xw;
      phi1   = -DEG2RAD*yh;
    }
    else
    {
      theta0= -atan((x-xw/2.0)/dist);              /* Angles to aim at target */
      phi0  = -atan((y-yh/2.0)/dist);
      theta1= -atan((x+xw/2.0)/dist);
      phi1  = -atan((y+yh/2.0)/dist);
    }

    /* shot towards target : flat distribution */

    if (gaussian)
    {
      theta= theta0+(theta1- theta0)*(randnorm()*FWHM2RMS+0.5);
      phi  = phi0  +(phi1  - phi0)  *(randnorm()*FWHM2RMS+0.5);
    }
    else
    {
      theta= theta0+(theta1- theta0)*rand01();
      phi  = phi0  +(phi1  - phi0)  *rand01();
    }
    /* Assume linear distribution */
    lambda = lambda_0+d_lambda*randpm1();
    if (lambda <= 0) ABSORB;

    v = K2V*(2*PI/lambda);

    p = p_in * fabs((theta1 - theta0)*(phi1 - phi0));

    if (flux_file && strlen(flux_file) > 0)
    {
       p *= exp(Table_Value(pTable, lambda, 1));
    }
    else if (T1 > 0 && I1 > 0)
    {
      lambda2 = lambda*lambda;
      lambda5 = lambda2*lambda2*lambda;
      Maxwell = I1 * L2P/lambda5 * exp(-lambda02/lambda2);  /* 1/AA */

      if ((T2 > 0) && (I2 > 0))
      {
        Maxwell += I2 * L2Pb/lambda5 * exp(-lambda02b/lambda2);
      }
      if ((T3 > 0) && (I3 > 0))
      {
        Maxwell += I3 * L2Pc/lambda5 * exp(-lambda02c/lambda2);
      }
      p *= Maxwell;
    }

    /* Perform the correct treatment - no small angle approx. here! */
    tan_h = tan(theta);
    tan_v = tan(phi);
    vz = v / sqrt(1 + tan_v*tan_v + tan_h*tan_h);
    vy = tan_v * vz;
    vx = tan_h * vz;
    SCATTER;
  }
%}

FINALLY
%{
  Table_Free(&pTable);
%}

MCDISPLAY
%{
  double xmin;
  double xmax;
  double ymin;
  double ymax;

  if ((h == 0) || (w == 0))
  {
    magnify("xy");
    circle("xy",0,0,0,radius);
    if (gaussian)
      circle("xy",0,0,0,radius/2);
  }
  else
  {
    xmin = -w/2; xmax = w/2;
    ymin = -h/2; ymax = h/2;

    magnify("xy");
    multiline(5, (double)xmin, (double)ymin, 0.0,
             (double)xmax, (double)ymin, 0.0,
             (double)xmax, (double)ymax, 0.0,
             (double)xmin, (double)ymax, 0.0,
             (double)xmin, (double)ymin, 0.0);
    if (gaussian)
      circle("xy",0,0,0,sqrt(w*w+h*h)/4);
  }
%}

END