#define VERSION 5
#define TYPE 4

#include <math.h>
#include "cplx.h"
#include <p_funct.h>

#define arrayMax 50                 /* The parameters' maximum size, adjust as needed */
#define PfileName "parameters.dat"  /* The parameter data file's name, adjust as needed */

double ErrFunc1(double y);
void initParam();
void calcReflec(double x, double *result);

float P[arrayMax], C[arrayMax];

setup() {
}

user4()
{
  double x, y,        /* The data x and y */
         xMax, xMin,  /* The x-axis ranges */
         interval;    /* An interval to increment x by */

  int i;              /* A counting number */

  Point pt;           /* A data point in C-Plot */

  initParam();        /* Initialize the parameters using the data file */

  xMax = get_xmax();
  xMin = get_xmin();
  interval = (xMax - xMin) / 200;

  printf("200 points will be generated and sent to C-Plot...\n");
  set_npts(200);

  for (x=xMin, i=0; x<=xMax; x+=interval, i++) {
    calcReflec(x, &y);
    set_x(&pt, x);
    set_y(&pt, y);
    pt_set(&pt, i);
  }

  printf("\nDone.\n");
}


/*****************************************************************************************
Error function definition from Abramowitz and Segun page 299 7.1.25 */

double ErrFunc1(double y)
{
  double yx, sgn, tx, zx;

  yx=sqrt(y*y)/sqrt(2.);
  sgn=1.0;

  if (y <= 0.0) sgn=-1.;

  tx=1./(1. + .47047*yx);
  zx= 1.-(.34802*tx-.0958798*tx*tx+.7478556*tx*tx*tx)*exp(-yx*yx); 
  zx =-0.5*sgn*zx;

  return zx;      
}

/****************************************************************************************/
void calcReflec(double x, double *result)
{
  struct cplx f[21], an[21], rn[21], fnmax, f1, rn1;
  double rho[21], d[21], my[21], sig, del[21], bet[21], xy[6], y;
  double slu, xsum, rsum, theta, frsnll, zx;
  int nmax, n, j;

  nmax      = C[1] + 2;
  sig       = P[2];
  rho[1]    = 0.0;
  my[1]     = 0.0;
  d[1]      = 0.0;
  rho[nmax] = P[4];
  my[nmax]  = P[5];
  d[nmax]   = 0.0;
  slu       = 1.0;  /* For neutrons. */

  if (C[7] == 0.0) slu = 2.817938e-5;     /* For X-rays. r-e in [AA] */

  if (C[2] == 0.0) {                            /* Refelectivity Calculations */
    theta = x*P[1]/(2.0*C[5]);

    del[nmax] = slu * 2.0 * 3.14159 / (C[5]*C[5]) * rho[nmax];
    bet[nmax] = .5/C[5]*1.e-8*my[nmax];
    f[nmax] = cplxsqrt(makecplx(theta*theta - 2.0 * del[nmax],-2.0 * bet[nmax]));

    del[1] = 0.0;
    bet[1] = 0.0;

    f[1] = makecplx(theta,0); 

    for (n=nmax;n>=3;n=n-1) {
      d[n-1]=P[3*(nmax-n)+6];         /*LAYER N-1 */
      rho[n-1]=P[3*(nmax-n)+7]; 
      my[n-1]=P[3*(nmax-n)+8];

      del[n-1]=slu*2.*3.14159/(C[5]*C[5])*rho[n-1];       
      bet[n-1]=.5/C[5]*1.e-8*my[n-1];
      f[n-1]=cplxsqrt(makecplx(theta*theta-2.0*del[n-1],-2.0*bet[n-1]));
    }

/* Calculate Fresnel law for subphase-only: */

    fnmax=cplxsqrt(makecplx(theta*theta-2.0*del[nmax],-2.0*bet[nmax]));
    f1=makecplx(theta,0.0);
    rn1= cdiv( cdiff(f1,fnmax),csum(f1,fnmax) );
    frsnll=modulus(rn1)*modulus(rn1);
    y=frsnll;

    rn[nmax]=makecplx(0.0,0.0);
       
    for (n=nmax; n>=2; n--) {
      an[n-1]=cdiv( cdiff (f[n-1],f[n]),csum(f[n-1],f[n]) );
      rn[n-1]=cmult( cplxexp( cmult(makecplx(0,-2.*C[5]*d[n-1]),f[n-1])), cdiv( csum(rn[n],an[n-1]),csum(cmult(rn[n],an[n-1]),makecplx(1.0,0) )));
    }

    if(C[8] > 0.0) {  /* Phase Output */ 
      *result=phase(rn[1]);
      return;
    }

/* Output Of Reflectivity: */

    y = P[3]*modulus(rn[1])*modulus(rn[1])*exp(- P[1]*P[1]*x*x*sig*sig);

    if (C[6] != 0.0) y = y/frsnll; 
    if (C[3] != 0.0) y = log10(y);

    *result=y;

    return;                                
  } /* End of Reflectivity Calculations */

  /* Real Space (Rho vs. X=Height: */

  if (C[4] > 0) {   /* SMAEARING */ 
    if (nmax > 2)
      for (n=nmax; n>=3; n--) {
	rho[n-1] = P[3 * (nmax-n) + 7];
	d[n-1] = P[3 * (nmax-n) + 6];
      }

    xy[1] = x/sig;
    for (j=2; j<=nmax-1; j++) {
      xy[j] = xy[j-1] - d[nmax + 1 - j]/sig;
    }

    rsum=0.0;
    
    for (j=nmax-1; j>=1; j--) {
      zx = xy[nmax-j];
      rsum = rsum + ErrFunc1(zx) * (rho[j+1]-rho[j]);
    }

    *result = rho[nmax] / 2 + rsum;

    return;
  }
}

/*********************************************************************************/
void initParam()
{
  FILE *inf;
  char ch, dummy;
  int n;       /* The array index */
  float val;   /* The array value */

  inf=fopen(PfileName, "r");     /* The location of the parameters */

  if (!inf) {
    printf("There was an error opening the file %s!\n", PfileName);
    exit(1);
  }

  do {
    ch=fgetc(inf);
    switch (toupper(ch)) {
    case '#':                           /* Ignore comment lines */
	do { dummy=fgetc(inf); } while (dummy!='\n' && dummy!=EOF);
	break;
    case 'P':                           /* This is the parameter P */
	while (fgetc(inf)!='[');        /* Skip all characters until a [ is reached */
	fscanf(inf, "%d", &n);          /* The array index */

	while (fgetc(inf)!='=');        /* Skip all characters until an = is reached */
	fscanf(inf, "%f", &val);        /* The array value */

	P[n]=val;
	break;
    case 'C':
	while (fgetc(inf)!='[');         /* Skip all characters until a [ is reached */
	fscanf(inf, "%d", &n);           /* The array index */

	while (fgetc(inf)!='=');         /* Skip all characters until an = is reached */
	fscanf(inf, "%f", &val);         /* The array value */

	C[n]=val;
	break;
    default:                             /* Ignore any other characters */
      break;
    }
  } while (ch!=EOF);       /* Continue until the end of file marker */

  fclose(inf);
}