/*******************************************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Source_extended
*
* %Identification
* Written by:Arne 'S Jegers
* Date: May 6, 2019
* Origin: Technical University of Denmark
* Release: McXtrace 1.4
*
* A plane source emitting x-rays to emulate a distant extended source, such as a
* nebula or galaxy
*
* %Description
* A rectangular x-ray source that samples ray intensity from an image, and
* deflects the emitted ray to reflect having been emitted from the sampled part
* of the extended source. Rays that are sampled from the same point on the image
* are collimated, but can be emitted from anywhere on the rectangular source.
* The image should be provided as a 2D-ascii table, whose header includes the
* following entries:
*
* w_pixels and h_pixels: The width and height of the image in pixels
* x_ref and y_ref: x and y reference values of the FITS image
* r_max: The maximum distance from the point (x_ref, y_ref) to any corner of the image
* iCD11, iCD12, iCD21 and iCD22: The values of the FITS image's CD matrix
*
* %Parameters
* yheight [m]   Height of rectangle in (x,y,0) plane where x-rays
*               are generated.
* xwidth  [m]   Width of rectangle in (x,y,0) plane where x-rays
*               are generated.
* E0:     [keV] Mean energy of xrays.
* dE:     [keV] Energy half spread of x-rays (flat or gaussian sigma).
* lambda0:[AA]  Mean wavelength of x-rays.
* dlambda:[AA]  Wavelength half spread of x-rays.
* gauss:  [1]   Gaussian (1) or Flat (0) energy/wavelength distribution
* flux: [pht/s] total flux radiated from the source
* incoherent: [] Source is fully incoherent
* phase: [] Set phase to something given.
* image_path: [string] Path to file containing the heat map of the source 
* %End
*******************************************************************************/

DEFINE COMPONENT Source_extended
SETTING PARAMETERS (string spectrum_file=NULL,yheight=0, xwidth=0,
  dist=0, E0=0, dE=0, lambda0=0, dlambda=0, flux=0,gauss=0,incoherent=1,phase=0, string image_path="")

SHARE
%{
  %include "read_table-lib"
  struct Sextended_prms {
    double l0, dl;
    double pmul, pint;
    t_Table T;
  };
  struct Sextended_table_prms {
    double xref, yref, rmax;
    int pw, ph;
    t_Table data;
    double iCD[2][2];
  };
  double*
  sphericalToCartesian (double r, double theta, double Phi, double* cartesian) {
    cartesian[0] = r * cos (theta) * cos (Phi);
    cartesian[1] = r * cos (theta) * sin (Phi);
    cartesian[2] = r * sin (theta);
  }

  double*
  cartesianToSpherical (double x, double y, double z, double* spherical) {
    spherical[0] = sqrt (x * x + y * y + z * z);
    spherical[1] = atan2 (z, sqrt (x * x + y * y));
    spherical[2] = atan2 (y, x);
  }
%}

DECLARE
%{
  double srcArea;
  int square;
  struct Sextended_prms prms;
  struct Sextended_table_prms extendedParams;
%}

INITIALIZE
%{
  square = 0;
  srcArea = 0;
  square = 1;
  srcArea = xwidth * yheight;

  if (srcArea <= 0) {
    printf ("Source_flat: %s: Source area is <= 0 !\n ERROR - Exiting\n", NAME_CURRENT_COMP);
    exit (0);
  }

  int status = 0;
  if (status = Table_Read (&(extendedParams.data), image_path, 0) == -1) {
    fprintf (stderr, "Source_extended(%s) Error: Could not parse file \"%s\"\n", NAME_CURRENT_COMP, image_path ? image_path : "");
    exit (-1);
  }

  char** header_parsed
      = Table_ParseHeader (extendedParams.data.header, "w_pixels=", "h_pixels=", "x_ref=", "y_ref=", "r_max=", "iCD11=", "iCD12=", "iCD21=", "iCD22=", NULL);
  if (header_parsed[0] && header_parsed[1] && header_parsed[2] && header_parsed[3] && header_parsed[4] && header_parsed[5] && header_parsed[6]
      && header_parsed[7]) {
    printf (header_parsed[0]);
    printf ("\n");
    extendedParams.pw = strtod (header_parsed[0], NULL);
    extendedParams.ph = strtod (header_parsed[1], NULL);
    extendedParams.xref = strtod (header_parsed[2], NULL);
    extendedParams.yref = strtod (header_parsed[3], NULL);
    extendedParams.rmax = strtod (header_parsed[4], NULL);
    extendedParams.iCD[0][0] = strtod (header_parsed[5], NULL);
    extendedParams.iCD[0][1] = strtod (header_parsed[6], NULL);
    extendedParams.iCD[1][0] = strtod (header_parsed[7], NULL);
    extendedParams.iCD[1][1] = strtod (header_parsed[8], NULL);
  }

  if (spectrum_file) {
    /*read spectrum from file*/
    int status = 0;
    if ((status = Table_Read (&(prms.T), spectrum_file, 0)) == -1) {
      fprintf (stderr, "Source_extended(%s) Error: 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;
    prms.pint = 0;
    t_Table* T = &(prms.T);
    for (i = 0; i < prms.T.rows - 1; i++) {
      prms.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]);
    }
    printf ("Source_flat(%s) Integrated intensity radiated is %g pht/s\n", NAME_CURRENT_COMP, prms.pint);
    if (E0)
      printf ("Source_flat(%s) E0!=0 -> assuming intensity spectrum is parametrized by energy [keV]\n", NAME_CURRENT_COMP);
  } else if (!E0 && !lambda0) {
    fprintf (stderr, "Source_flat(%s): Error: Must specify either wavelength or energy distribution\n", NAME_CURRENT_COMP);
    exit (0);
  }
  /*calculate the X-ray weight from the flux*/
  if (flux) {
    prms.pmul = flux;
  } else {
    prms.pmul = 1;
  }
  prms.pmul *= 1.0 / ((double)mcget_ncount ());
%}


TRACE
%{
  double chi, e, k, l, r, pdir;
  char done = 0;
  double xpix, ypix, theta, Phi;
  double xDeg, yDeg, xRel, yRel;

  Phi = 0;
  z = 0;
  if (square == 1) {
    x = xwidth * (rand01 () - 0.5);
    y = yheight * (rand01 () - 0.5);
  }

  /*pdir contains the unnormalized solid angle weighting */
  p = 1;

  if (spectrum_file) {
    double pp = 0;
    // while (pp<=0){
    l = prms.T.data[0] + (prms.T.data[(prms.T.rows - 1) * prms.T.columns] - prms.T.data[0]) * rand01 ();
    pp = Table_Value (prms.T, l, 1);
    //}
    p *= pp;
    /*if E0!=0 convert the tabled value to wavelength*/
  } else if (E0) {
    if (!dE) {
      e = E0;
    } else if (gauss) {
      e = E0 + dE * randnorm ();
    } else {
      e = randpm1 () * dE * 0.5 + 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);
  }

  // ===========================================================================
  //  This section randomly picks a direction to emit a ray to, and samples the
  //  source image for an intensity at the appropriate point

  while (!done) {
    theta = sqrt (rand01 ()) * extendedParams.rmax / 180 * M_PI; // normalized
    Phi = rand01 () * 2 * M_PI;

    xDeg = tan (theta) * cos (Phi) * 180 / M_PI;
    yDeg = tan (theta) * sin (Phi) * 180 / M_PI;

    xRel = extendedParams.iCD[0][0] * xDeg + extendedParams.iCD[0][1] * yDeg;
    yRel = extendedParams.iCD[1][0] * xDeg + extendedParams.iCD[1][1] * yDeg;

    xpix = xRel + extendedParams.xref;
    ypix = yRel + extendedParams.yref;

    if (xpix > 0 && xpix < extendedParams.pw && ypix > 0 && ypix < extendedParams.ph) {
      done = 1;
    }
  }
  // Both are inverted because, relative to the spherical coordinate system of
  // theta and Phi, the rays are emitted 'backwards'
  kx = -sin (theta) * cos (Phi) * k;
  ky = -sin (theta) * sin (Phi) * k;

  kz = cos (theta) * k;

  p = Table_Value2d (extendedParams.data, ypix, xpix);

  /*randomly pick phase or set to something real*/
  if (incoherent) {
    Phi = rand01 () * 2 * M_PI;
  } else {
    Phi = phase;
  }

  /*set polarization vector*/
  Ex = 0;
  Ey = 0;
  Ez = 0;
%}

FINALLY
%{
  Table_Free (&(prms.T));
  Table_Free (&(extendedParams.data));
%}

MCDISPLAY
%{
  if (square == 1) {
    magnify ("xy");
    rectangle ("xy", 0, 0, 0, xwidth, yheight);
  }
%}

END