/*******************************************************************************
*
* 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