// Parratt approximation of a structure factor
// parameters A1-A7 used to define a density function, which you can modify as you like


#include <origin.h>
#include <math.h>

#include "rho.cpp"
//#include "fresnel.cpp"

// parameters below define number of slices and starting and finishing point of the density
// set it to cover whole interval in which your density function changes, and try 
// different PMAXINT values to be sure the function converges
double PLAMBDA=1.534;
double f0f1z=1;
double f2z=0;

double PMIN=-20;
double PMAX=100;
double PMAXINT=500;

double parratt(double Qz,double A1,double A2,double A3,double A4,double A5,double A6,double A7,double A8)
{
	
complex n_j[4000];
complex k_j[4000];
complex r_j[4000];
complex x_j;
complex undere1,undere2;

double delta,beta,fr,rr;
double step,Qc;
double k;
double result;

double L22P;
double rho_e=2.814e-5;

k=2.0*pi/PLAMBDA;
L22P=rho_e*PLAMBDA*PLAMBDA/(2.0*pi);
step=(PMAX-PMIN)/PMAXINT;
Qc=4.0*sqrt(pi*rho_e*rho(PMAX,A1,A2,A3,A4,A5,A6,A7,A8));

for(int ii=0;ii<=PMAXINT;ii++)
	{	rr=rho(PMIN+step*ii,A1,A2,A3,A4,A5,A6,A7,A8);
		n_j[ii].m_re=1.0-L22P*f0f1z*rr;
		n_j[ii].m_im=L22P*f2z*rr;
		k_j[ii]=0.5*sqrt(4.0*k*k*(n_j[ii]*n_j[ii]-1.0)+Qz*Qz);
	}

for(ii=0;ii<PMAXINT;ii++)
	r_j[ii]=(k_j[ii]-k_j[ii+1])/(k_j[ii]+k_j[ii+1]);

x_j.m_re=0.0;
x_j.m_im=0.0;
for(ii=PMAXINT-1;ii>=0;ii--)
	{   
		undere1.m_re=0.0;
		undere1.m_im=2.0*(PMIN+step*ii);
		undere2=-k_j[ii]*undere1;
		undere1*=k_j[ii+1];
		x_j=exp(undere2)*((x_j*exp(undere1)+r_j[ii])/(x_j*exp(undere1)*r_j[ii]+1.0));
	}
result=cabs(x_j);
result*=result;
if(Qc<Qz)
	fr=((Qz-sqrt(Qz^2-Qc^2))/(Qz+sqrt(Qz^2-Qc^2)))^2;
else
	fr=((Qz-sqrt(Qc^2-Qz^2))/(Qz+sqrt(Qc^2-Qz^2)))^2;
result/=fr;//Rf(Qz,Qc,0,PLAMBDA);
return(result);	
}


// parameters below define "integration" interval and number of Z-intervals, 
// set it to cover whole interval in which your density function changes, and try 
// different MAXINT values to be sure the function converges
double ZMAX=100;
double ZMIN=-20;
double MAXINT=500;

double born(double Qz,double A1,double A2,double A3,double A4,double A5,double A6,double A7,double A8)
{
complex sum,xx;
double res,step;

sum.m_re=0.0;
sum.m_im=0.0;
step=(ZMAX-ZMIN)/MAXINT;
for(int ii=0;ii<MAXINT-1;ii++)
	{	
		xx.m_re=0.0;
		xx.m_im=Qz*(ZMIN+step*ii);
		xx=exp(xx);
		sum+=xx*(rho(ZMIN+step*(ii+1),A1,A2,A3,A4,A5,A6,A7,A8)-rho(ZMIN+step*(ii),A1,A2,A3,A4,A5,A6,A7,A8));
	}
sum/=rho(ZMAX,A1,A2,A3,A4,A5,A6,A7,A8);
res=cabs(sum);
res*=res;
return(res);
}