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

int IDENSITY=2; 
// 0 - user defined
// 1 - distored crystal model
// 2 - two layers on a substrate


//DENSITY FUNCTION
double rho(double zz,double A1,double A2,double A3,double A4,double A5,double A6,double A7,double A8)
{


double result;
	
// USER DEFINED FUNCTION
if(IDENSITY==0)
{
// You can type in your own function using A1-A7 parameters
result=1.0;

return(result);
}


// STANDART DISTORED CRYSTAL MODEL
// A1 - number of summation steps
// A2 - sigma average
// A3 - sigma0
// A4 - spacing
// A5 - z0 offset
// A6 - rho in the bulk
// A7 - not used
if(IDENSITY==1)
{
double sn;

result=0.0;
for(int ii=0;ii<A1;ii++)
	{	
		sn=sqrt(ii*A2*A2+A3*A3);
		result+=(A4/(sn*sqrt(2.0*pi)))*exp(-(zz-ii*A4-A5)^2/(2.0*sn*sn));
	}
return(A6*result);
}

// TWO LAYERS ON SUBSTRATE
// A1 - 1st layer thickness
// A2 - surface rhoughness
// A3 - 1 st layer density
// A4 - 2nd layer thickness
// A5 - 1st-to-2nd layer roughtness
// A6 - 2nd layer density
// A7 - substrate roughness
// A8 - substrate density
if(IDENSITY==2)
{
result=0.0;
double s2=1.41421356;

result=A3*(erf(zz/(s2*A2))+1.0)/2.0+A6*(erf((zz-A1)/(s2*A5))+1.0)/2.0-A3*(erf((zz-A1)/(s2*A5))+1.0)/2.0+\
A8*(erf((zz-A1-A4)/(s2*A7))+1.0)/2.0-A6*(erf((zz-A1-A4)/(s2*A7))+1.0)/2.0;
return(result);
}

if(IDENSITY==3)
{
result=0.0;

result=0.6679*(erf(zz/3.03)+1.0)/2.0+2.104*(erf((zz-6.45)/3.5637)+1.0)/2.0-0.6679*(erf((zz-6.45)/3.5637)+1.0)/2.0\
+0.52236*(erf((zz-18.1612-6.45)/5.3773)+1.0)/2.0-2.104*(erf((zz-18.1612-6.45)/5.3773)+1.0)/2.0+\
1*(erf((zz-16.7353-18.1612-6.45)/6.5329)+1.0)/2.0-0.52236*(erf((zz-16.7353-18.1612-6.45)/6.5329)+1.0)/2.0;;
return(result);
}
}