/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Source_gen
*
* %I
* Written by: Emmanuel Farhi, Kim Lefmann, edits Jonas O Birk
* Date: Aug 27, 2001
* Version: $Revision: 1.41 $
* Origin: ILL/Risoe
* Release: McStas 1.12c
* 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 or spatially gaussian neutron source with flat or 
* Maxwellian energy/wavelength spectrum.
*
* %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). A high-energy tail of shape 
* A/(lambda0-lambda) is added by setting the HEtailA and HEtailL0 parameters.
* 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.
* Additionl distributions for the horizontal and vertical phase spaces
* distributions (position-divergence) may be specified with the
* 'xdiv_file' and 'ydiv_file' parameters. Format is free text, requiring
* a comment line '# xylimits: pos_min pos_max div_min div_max' to set
* the axis of the distribution matrix. All these files may be generated using
* standard monitors (better in McStas/PGPLOT format), e.g.:
*   Monitor_nD(options="auto lambda per cm2")
*   Monitor_nD(options="x hdiv, all auto")
*   Monitor_nD(options="y vdiv, all auto")
* 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 divergence uniform,.
* 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,I1=1e13,xw=0.01,yh=0.01)
*   Source_gen(h=0.1,w=0.1,Emin=1,Emax=3,I1=1e13,verbose=1,xw=0.01,yh=0.01)
*   EXTEND
*   %{
*      t = rand0max(1e-3); // set time from 0 to 1 ms for TOF instruments.
*   %}
*
* <b>Some sources parameters:</b>
* PSI cold source T1= 296.16, I1=8.5E11, T2=40.68, I2=5.2E11
* 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 = full horz. div  in deg
* yh:     [m]   Height(y) of target.
*                 If dist=0.0, yh = full 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
* HEtailA:  [1/(cm**2*st*AA**2)] Amplitude of heigh energy tail
* HEtailL0  [AA]  Offset of heigh energy tail
flux_file:[str] Name of a two columns [lambda flux] text file that contains
*                 the wavelength distribution of the flux in [1/(s*cm**2*st*AA)].
*                 Comments (#) and further columns are ignored. Format is
*                 compatible with McStas/PGPLOT wavelength monitor files.
*                 When specified, temperature and intensity values are ignored.
* flux_file_perAA:[1] When true (1), indicates that flux file data is already per
*                 Angstroem. If false, file data is per wavelength bin.
* flux_file_log:[1]   When true, will transform the flux table in log scale to
*                 improve the sampling. This may also cause
* xdiv_file:[str] Name of the x-horiz. divergence distribution file, given as a
*                 free format text matrix, preceeded with a line
*                   '# xylimits: xmin xmax xdiv_min xdiv_max'
* ydiv_file:[str] Name of the y-vert. divergence distribution file, given as a
*                 free format text matrix, preceeded with a line
*                   '# xylimits: ymin ymax ydiv_min ydiv_max'
* verbose:  [0/1] display info about the source. -1 unactivate source.
*
* %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 (string flux_file=0, string xdiv_file=0, string ydiv_file=0)
SETTING PARAMETERS (radius=0.0, dist=0, xw=0, yh=0, E0=0, dE=0, Lambda0=0, dLambda=0, I1=0,
                    h=0, w=0, verbose=0, T1=0,
                    flux_file_perAA=0, flux_file_log=0,
                    Lmin=0,Lmax=0,Emin=0,Emax=0,T2=0,I2=0,T3=0,I3=0,length=0,
		    HEtailA=0,HEtailL0=0)
OUTPUT PARAMETERS (p_in, lambda0, lambda02, L2P, lambda0b, lambda02b, L2Pb,lambda0c, lambda02c, L2Pc, pTable, pTable_x, pTable_y,pTable_xmin, pTable_xmax, pTable_xsum, pTable_ymin, pTable_ymax, pTable_ysum, pTable_dxmin, pTable_dxmax, pTable_dymin, pTable_dymax)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

SHARE
%{
%include "read_table-lib"

#ifndef SOURCE_GEN_DEF
#define SOURCE_GEN_DEF
/*******************************************************************************
* str_dup_numeric: makes a clean copy of a string and allocate as numeric
*******************************************************************************/
char *str_dup_numeric(char *orig)
{
  long i;
  char *valid;

  if (!orig || !strlen(orig)) return(NULL);

  for (i=0; i < strlen(orig); i++)
  {
    if ( (orig[i] > 122)
      || (orig[i] < 32)
      || (strchr("!\"#$%&'()*,:;<=>?@[\\]^`/ ", orig[i]) != NULL) )
    {
      orig[i] = ' ';
    }
  }
  orig[i] = '\0';
  /* now skip spaces */
  for (i=0; i < strlen(orig); i++) {
    if (*orig == ' ') orig++;
    else break;
  }

  return(orig);
} /* str_dup_numeric */
#endif
%}

DECLARE
%{
 double p_in;
 double Lambda0=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;
 t_Table pTable_x;
 t_Table pTable_y;
 double  pTable_xmin;
 double  pTable_xmax;
 double  pTable_xsum=0;
 double  pTable_ymin;
 double  pTable_ymax;
 double  pTable_ysum=0;
 double  pTable_dxmin;
 double  pTable_dxmax;
 double  pTable_dymin;
 double  pTable_dymax;

%}

INITIALIZE
%{
  double source_area, k;

  /* spectrum characteristics */
  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));
    /* put table in Log scale */
    int i;
    if (pTable.columns < 2) exit(fprintf(stderr, "Source_gen: %s: Flux file %s should contain at least 2 columns\n", NAME_CURRENT_COMP, flux_file));
    double table_lmin=FLT_MAX, table_lmax=-FLT_MAX;
    double tmin=FLT_MAX, tmax=-FLT_MAX;
    for (i=0; i<pTable.rows; i++) {
      double val = Table_Index(pTable, i,1);
      val = Table_Index(pTable, i,0); /* lambda */
      if (val > tmax) tmax=val;
      if (val < tmin) tmin=val;
    }
    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);
      if (flux_file_log)
        val = log(val > 0 ? val : tmin/10);
      Table_SetElement(&pTable, i, 1, val);
      val = Table_Index(pTable, i,0); /* lambda */
      if (val > table_lmax) table_lmax=val;
      if (val < table_lmin) table_lmin=val;
    }
    if (!Lmin && !Lmax && !Lambda0 && !dLambda && !E0 && !dE && !Emin && !Emax) {
      Lmin = table_lmin; Lmax = table_lmax;
    }
    if (Lmax > table_lmax) {
      if (verbose) fprintf(stderr, "Source_gen: %s: Maximum wavelength %g is beyond table range upper limit %g. Constraining.\n", NAME_CURRENT_COMP, Lmax, table_lmax);
      Lmax = table_lmax;
    }
    if (Lmin < table_lmin) {
      if (verbose) fprintf(stderr, "Source_gen: %s: Minimum wavelength %g is below table range lower limit %g. Constraining.\n", NAME_CURRENT_COMP, Lmin, table_lmin);
      Lmin = table_lmin;
    }
  } 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 = Lambda0; }

    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 = Lambda0; }

    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 = Lambda0; }
  }

  /* now read position-divergence files, if any */
  if (xdiv_file && strlen(xdiv_file) > 0) {
    int i,j;
    if (Table_Read(&pTable_x, xdiv_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, xdiv_file));
    pTable_xsum = 0;
    for (i=0; i<pTable_x.rows; i++)
      for (j=0; j<pTable_x.columns; j++)
        pTable_xsum += Table_Index(pTable_x, i,j);

    /* now extract limits */
    char **parsing;
    char xylimits[1024];
    strcpy(xylimits, "");
    parsing = Table_ParseHeader(pTable_x.header,
      "xlimits", "xylimits",
      NULL);

    if (parsing) {
      if (parsing[0])  strcpy(xylimits, str_dup_numeric(parsing[0]));
      if (parsing[1] && !strlen(xylimits))
                       strcpy(xylimits, str_dup_numeric(parsing[1]));
      for (i=0; i<=1; i++) {
        if (parsing[i]) free(parsing[i]);
      }
      free(parsing);
    }
    i = sscanf(xylimits, "%lg %lg %lg %lg",
      &(pTable_xmin),  &(pTable_xmax),
      &(pTable_dxmin), &(pTable_dxmax));
    if (i != 2 && i != 4 && verbose)
      fprintf(stderr, "Source_gen: %s: invalid xylimits '%s' from file %s. extracted %i values\n",
        NAME_CURRENT_COMP, xylimits, xdiv_file, i);

    if (!w) w=pTable_xmax-pTable_xmin;
    if (!xw && !dist) xw=fabs(pTable_dxmax-pTable_dxmin);
  }

  if (ydiv_file && strlen(ydiv_file) > 0) {
    int i,j;
    if (Table_Read(&pTable_y, ydiv_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, ydiv_file));
    pTable_ysum = 0;
    for (i=0; i<pTable_y.rows; i++)
      for (j=0; j<pTable_y.columns; j++)
        pTable_ysum += Table_Index(pTable_y, i,j);

    /* now extract limits */
    char **parsing;
    char xylimits[1024];
    strcpy(xylimits, "");
    parsing = Table_ParseHeader(pTable_y.header,
      "xlimits", "xylimits",
      NULL);

    if (parsing) {
      if (parsing[0])  strcpy(xylimits,str_dup_numeric(parsing[0]));
      if (parsing[1] && !strlen(xylimits))
                       strcpy(xylimits,str_dup_numeric(parsing[1]));
      for (i=0; i<=1; i++) {
        if (parsing[i]) free(parsing[i]);
      }
      free(parsing);
    }
    i = sscanf(xylimits, "%lg %lg %lg %lg",
      &(pTable_ymin),  &(pTable_ymax),
      &(pTable_dymin), &(pTable_dymax));
    if (i != 2 && i != 4 && verbose)
      fprintf(stderr, "Source_gen: %s: invalid xylimits '%s' from file %s. extracted %i values\n",
        NAME_CURRENT_COMP, xylimits, ydiv_file, i);
    if (!h)  h=pTable_ymax-pTable_ymin;
    if (!yh && !dist) yh=fabs(pTable_dymax-pTable_dymin);
  }

  /* tests for parameter values */
  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))
  { 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);
  Lambda0=fabs(Lambda0); dLambda=fabs(dLambda);
  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 full divergence in [deg]\n",
                   NAME_CURRENT_COMP);
  }
  Lmin = Lambda0 - dLambda; /* Angstroem */
  Lmax = Lambda0 + dLambda;

  /* compute initial weight factor p_in to get [n/s] */
  if ((I1 > 0  && T1 >= 0) || (flux_file && strlen(flux_file) > 0))
  { /* the I1,2,3 are usually in [n/s/cm2/st/AA] */
    if (radius)
      source_area = radius*radius*PI*1e4; /* circular cm^2 */
    else
      source_area = h*w*1e4; /* square cm^2 */
    p_in  = source_area; /* cm2 */
    p_in *= (Lmax-Lmin); /* AA. 1 bin=AA/n */
    if (flux_file && strlen(flux_file) && !flux_file_perAA)  p_in *= pTable.rows/(Lmax-Lmin);
  }
  else
    p_in = (I1 > 0? I1 : 1)/4/PI; /* Small angle approx. */
  p_in /= mcget_ncount();
  if (!T1 && I1) p_in *= I1;

  if (radius == 0 && h == 0 && w == 0)
  {
    fprintf(stderr,"Source_gen: %s: Error: Please specify source geometry (radius, h, w)\n",
      NAME_CURRENT_COMP);
    exit(-1);
  }
  if (xw*yh == 0)
  {
    fprintf(stderr,"Source_gen: %s: Error: Please specify source target (xw, yh)\n",
      NAME_CURRENT_COMP);
    exit(-1);
  }

  if (verbose)
  {
    printf("Source_gen: component %s ", NAME_CURRENT_COMP);
    if ((h == 0) || (w == 0))
      printf("(disk, radius=%g)", radius);
    else
      printf("(square %g x %g)",h,w);
    printf("\n            spectra ");
    printf("%.3f to %.3f AA (%.3f to %.3f meV)", Lmin, Lmax, 81.81/Lmax/Lmax, 81.81/Lmin/Lmin);
    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>=0 && 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/cm2].\n");
    }
    else
    { printf("  Flux is Phi in [n/s].\n");
    }
    if (xdiv_file && strlen(xdiv_file) > 0)
      printf("  File %s x=[%g:%g] [m] xdiv=[%g:%g] [deg] used as horizontal phase space distribution.\n", xdiv_file, pTable_xmin, pTable_xmax, pTable_dxmin, pTable_dxmax);
    if (ydiv_file && strlen(ydiv_file) > 0)
      printf("  File %s y=[%g:%g] [m] ydiv=[%g:%g] [deg] used as vertical phase space distribution.\n", ydiv_file, pTable_ymin, pTable_ymax, pTable_dymin, pTable_dymax);
  }
  else
    if (verbose == -1)
      printf("Source_gen: component %s unactivated", NAME_CURRENT_COMP);
%}

TRACE
%{
  double dx,dy,xf,yf,rf,pdir,chi,v,r, lambda;
  double Maxwell, lambda2, lambda5;
  double xwg, yhg;

  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;
  /* Assume linear wavelength distribution */
    lambda = Lambda0+dLambda*randpm1();
    if (lambda <= 0) ABSORB;
    v = K2V*(2*PI/lambda);
      
    if (dist>0) {
      randvec_target_rect_real(&xf, &yf, &rf, &pdir,
			       0, 0, dist, xw, yh, ROT_A_CURRENT_COMP, x, y, z, 2);
      

      dx = xf-x;
      dy = yf-y;
      rf = sqrt(dx*dx+dy*dy+dist*dist);

      vz=v*dist/rf;
      vy=v*dy/rf;
      vx=v*dx/rf;
    } else {
      randvec_target_rect_angular(&vx, &vy, &vz, &pdir,
          0, 0, 1, xw*DEG2RAD, yh*DEG2RAD, ROT_A_CURRENT_COMP);
      vx *= v; vy *= v; vz *= v;
    }
    p = p_in*pdir;
    
    /* spectral dependency from files or Maxwellians */
    if (flux_file && strlen(flux_file) > 0)
    {
       double w=Table_Value(pTable, lambda, 1);
       if (flux_file_log) w=exp(w);
       p *= w;
    }
    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);
      }
      if (HEtailA>0)
      {
        Maxwell+=HEtailA/(lambda-HEtailL0)/(lambda-HEtailL0);
      }
      p *= Maxwell;
    }

    /* optional x-xdiv and y-ydiv weightening: position=along columns, div=along rows */
    if (xdiv_file && strlen(xdiv_file) > 0 && pTable_xsum > 0) {
      double i,j;
      j = (x-            pTable_xmin) /(pTable_xmax -pTable_xmin) *pTable_x.columns;
      i = (atan2(dx,rf)*RAD2DEG-pTable_dxmin)/(pTable_dxmax-pTable_dxmin)*pTable_x.rows;
      r = Table_Value2d(pTable_x, i,j); /* row, column */
      p *= r/pTable_xsum;
    }
    if (ydiv_file && strlen(ydiv_file) > 0 && pTable_ysum > 0) {
      double i,j;
      j = (y-            pTable_ymin) /(pTable_ymax -pTable_ymin) *pTable_y.columns;
      i = (atan2(dy,rf)*RAD2DEG-  pTable_dymin)/(pTable_dymax-pTable_dymin)*pTable_y.rows;
      r = Table_Value2d(pTable_y, i,j);
      p *= r/pTable_ysum;
    }
    SCATTER;
  }
%}

FINALLY
%{
  Table_Free(&pTable);
  Table_Free(&pTable_x);
  Table_Free(&pTable_y);
%}

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

  if ((h == 0) || (w == 0))
  {
    magnify("xy");
    circle("xy",0,0,0,radius);
  }
  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);
  }
%}

END