/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: curved_mirror
*
* %I
* Written by: Henrik Jacobsen
* Date: April 2012
* Origin: RNBI
*
* Single mirror plate that is curved and fits into an elliptic guide. 
*
* %D
*
*
* %P
* INPUT PARAMETERS:
*
* focus_s: [m]                 first focal point of ellipse in ABSOLUTE COORDINATES
* focus_e: [m]                 second focal point of ellipse in ABSOLUTE COORDINATES
* mirror_start: [m]            start of mirror in ABSOLUTE COORDINATES
* guide_start: [m]             start of guide in ABSOLUTE COORDINATES
* yheight: [m]                 the height of the mirror when not in the guide
* smallaxis: [m]               the smallaxis of the guide
* length: [m]                  the length of the mirror
* R0=0.99:		[]	low angle reflectivity
* Qc: [AA-1]                   critical scattering vector
* alpha: [AA]                  slope of reflectivity
* m:			[]	m-value of material
* W: [AA-1]                    width of supermirror cutoff
* transmit: [0/1]              0: non reflected neutrons are absorbed. 1: non reflected neutrons pass through
* substrate_thickness: [m]     thickness of substrate (for absorption)
* coating_thickness: [m]       thickness of coating (for absorption)
* theta_1:		[deg]	angle of mirror wrt propagation direction at start of mirror
* theta_2:		[deg]	angle of mirror wrt propagation direction at center of mirror
* theta_3:		[deg]	angle of mirror wrt propagation direction at end of mirror
* reflect: [q(Angs-1) R(0-1)]  (str)  Name of relfectivity file. Format 
*
* %D
*
* %E
*******************************************************************************/

DEFINE COMPONENT Mirror_Curved_Bispectral

SETTING PARAMETERS (string reflect=0,focus_s,focus_e, mirror_start, guide_start, yheight, smallaxis, length ,m, transmit=0, substrate_thickness=0.0005, coating_thickness=10e-6, theta_1=1.25, theta_2=1.25, theta_3=1.25)


SHARE
%{
%include "read_table-lib"
%}

DECLARE
%{
t_Table pTable;

%}

INITIALIZE
%{
  if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) {
    if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */
      exit(fprintf(stderr,"Mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect));
  }
%}

TRACE
%{
double f; //half of distance between focal points
double asquared;
double a; //half of ellipse length
double b; //half of ellipse width

double xprime; //x in coordinates with center of ellipse at (xprime,zprime)=(0,0)
double ymirror; //height of the mirror


//Defining the mirror
double a1;
double b1;
double c1;

//solving the time the neutron hits the sample
double A, B, C, D, E, P, Q, R, U, V, I, J, K;

//finding rotation of mirror
double alpha1, beta1, gamma1;
double theta_m;
double sin_theta_m, cos_theta_m;

double tan_theta_1;
double tan_theta_2;
double tan_theta_3;


double v_n; //speed of neutron perpendicular to surface

double Ref; //reflectivity

double dt;
double q;
  int intersect;

double discriminant;




double dt_2;
double dt_3;
int prop_case;
double x_2;
double y_2;
double z_2;
double t_2;
double x_3;
double y_3;
double z_3;
double t_3;

int x_hit;
int x_hit_2;
int x_hit_3;
double xprime_2;
double ymirror_2;
double xprime_3;
double ymirror_3;
int intersect_2;
int intersect_3;


intersect=0;
x_hit=0;
x_hit_2=0;
x_hit_3=0;
intersect_2=0;
intersect_3=0;
prop_case=0;

//printf("\n\n\n");
   double old_x = x, old_y = y, old_z = z, old_t=t, old_vx=vx, old_vz=vz, old_vy=vy;

// printf("x=%f, y=%f, z=%f, vx=%f, vy=%f, vz=%f\n",x,y,z,vx,vy,vz);

// Check if neutron hits mirror. First find which z,x coordinates it hits.

//mirror is defined by z(x)=a1x^3+b1x^2+c1x+d1, with dz/dx|x=-length/2=tan(theta_1), dz/dx|x=0=tan(theta_2), dz/dx|x=length/2=tan(theta3), z(0)=0. (d1=0)

tan_theta_1=tan(theta_1*DEG2RAD);
tan_theta_2=tan(theta_2*DEG2RAD);
tan_theta_3=tan(theta_3*DEG2RAD);


a1=2.0/3.0*(tan_theta_1+tan_theta_3-2.0*tan_theta_2)/(length*length);
b1=(tan_theta_3-tan_theta_1)/(2.0*length);
c1=tan_theta_2;


//neutron trajectory is defined by x=x0+vx*t, z=z0+vz*t. setting z=a1*x^3+b1*x^2+c1*x gives the equation A*t^3+B*t^2+C*t+D=0, with
A=a1*vx*vx*vx;
B=3.0*a1*x*vx*vx+b1*vx*vx;
C=3.0*a1*x*x*vx+2.0*b1*x*vx+c1*vx-vz;
D=a1*x*x*x+b1*x*x+c1*x-z;

//printf("a1=%f,b1=%f,c1=%f",a1,b1,c1);

//this equation must now be solved for t;

if (A!=0){
P=1/3.0*(3.0*C/A-B*B/(A*A));
Q=1/27.0*(2.0*B*B*B/(A*A*A)-9.0*B*C/(A*A)+27.0*D/A);

E=P*P*P/27.0+Q*Q/4.0;

// printf("A=%f, B=%f, C=%f, D=%f, 1e6P=%f, 1e6Q=%f, 1e6E=%f\n", A, B, C, D, 1e6*P, 1e6*Q, 1e6*E);

prop_case=0;
if (E>=0){

U=cbrt(-Q/2.0+sqrt(E));
V=cbrt(-Q/2.0-sqrt(E));

I=U+V-B/(3.0*A);
dt=I;
dt_2=I;
dt_3=I;
// printf("I=%f\n",I);

// J=-(U+V)/2+1i*(U-V)*sqrt(3)/2-B/(3*A) //complex solution
// K=-(U+V)/2-1i*(U-V)*sqrt(3)/2-B/(3*A) //complex solution
}else{
    R=acos(-Q/(2.0*sqrt(-P*P*P/27.0)));

// printf("R=%f\n",R);

   
   I=2.0*sqrt(fabs(P)/3.0)*cos(R/3.0)-B/A/3.0;
   J=-2.0*sqrt(fabs(P)/3.0)*cos(R/3.0 + 3.1415926535/3.0)-B/A/3.0;
   K=-2.0*sqrt(fabs(P)/3.0)*cos(R/3.0 - 3.1415926535/3.0)-B/A/3.0;

// printf("2.0*sqrt(abs(P)/3.0)=%f", 2.0*sqrt(abs(P)/3.0));
// printf("cos(R/3.0)=%f, cos(R/3.0 + 3.1415926535/3.0)=%f, cos(R/3.0 - 3.1415926535/3.0)=%f, -B/A/3.0=%f\n", cos(R/3.0), cos(R/3.0 + 3.1415926535/3.0), cos(R/3.0 - 3.1415926535/3.0), -B/A/3.0);

// printf("I=%f, J=%f, K=%f, \n",I, J, K);
// printf("P=%f, R=%f, A=%f, B=%f, \n",P, R, A, B);



// Three solutions. Find the smallest positive of these.
//there are problems with the solutions....
	if (I<=0){
		if (J<=0 && K<=0){dt=-1.0;} //if all three are negative, dt<0 and nothing happens
		if (J<=0 && K>0){dt=K;}  //if only K>0, dt=K
		if (J>0 && K<=0){dt=J;} //if only J>0, dt=J

		if (J>0 && K>0){	//if both J>0 and K>0, compare
			if (J>=K){dt=K; prop_case=1; dt_2=J;}else{dt=J; dt_2=I; prop_case=2;} } //dt is the smallest value
	}else{ //end if (I<=0)
		if (J<=0 && K<=0){dt=I;} //if only I>0, dt=I;

		if (J<=0 && K>0){ //if both I>0 and K>0, compare
			if (K>=I){dt=I; dt_2=K; prop_case=3;}else{dt=K; dt_2=I; prop_case=4;} } //dt is the smallest value

		if (J>0 && K<=0){ //if both I>0 and J>0, compare
			if (J>=I){dt=I; dt_2=J; prop_case=5;}else{dt=J; dt_2=I; prop_case=6;} } //dt is the smallest value

		if (J>0 && K>0){ //if all three>0, compare
			if (J>=K){ //either K or I is smallest
				if (K>=I){dt=I; if(J>=K){ dt_2=K; dt_3=J; prop_case=9;}else{dt_2=K; dt_3=J; prop_case=15;}}else{dt=K; if (J>=I){dt_2=I; dt_3=J; prop_case=10;}else{dt_2=J; dt_3=I; prop_case=11;} } //if K is smallest, compare it to I  
			}else{
				if (J>=I){dt=I; if (K>J){dt_2=J; dt_3=K; prop_case=12;}else{dt_2=J; dt_3=J; prop_case=16;}}else{dt=J; if (K>I){dt_2=I; dt_3=K; prop_case=13;}else{{dt_2=K; dt_3=I; prop_case=14;}}  }}  //else compare J to I
			} //end if(J>0 && K>0)
				
	} //end }else{ for if(I<=0)



}    // end }else{ for if (E>=0)
  

}else{ //end if (A!=0)  
if (B!=0){

discriminant=C*C-4*B*D;

if (discriminant<0){dt=-1.0;}else{ //only complex solutions: set dt<0 to avoid interaction
I=(-C-sqrt(discriminant))/(2.0*B);
J=(-C+sqrt(discriminant))/(2.0*B);

if (I<=0 && J<=0){dt = -1.0;} //both times are negative.
if (I<=0 && J>0 ){dt = J;} //set dt to only positive value.
if (I>0  && J<=0){dt = I;} //set dt to only positive value.
if (I>0  && J>0 ){if (I>J) {dt=J; dt_2=I; prop_case=7;}else{dt=I; dt_2=J; prop_case=8;} } //set dt to smallest positive value  

} //end if (discriminant<0){}else{
}else{ //end if (B!)=0
if (C!=0) { dt = -D/C;}else{
 printf("warning: A=B=C=0. Neutron is ignored\n"); }
} //end if(B!=0){}else{
} //end if (A!=0){}else{
//now intersection time has been found.

if (dt>0) { //if time is positive, propagate neutron to where it hits mirror. This is done without gravity.
// printf("before anything: x=%f,y=%f,z=%f,vx=%f,vy=%f,vz=%f, dt=%f\n",x,y,z,vx,vy,vz,dt);

    x += vx*dt;
    y += vy*dt;
    z += vz*dt;
    t += dt;


x_hit=(x >=-length/2 && x<=length/2);


if (prop_case==0){
x_2=x;
y_2=y;
z_2=z;
t_2=t;
x_3=x;
y_3=y;
z_3=z;
t_3=t;
}

if (prop_case>0)
{
x_2=old_x+vx*dt_2;
y_2=old_y+vy*dt_2;
z_2=old_z+vz*dt_2;
t_2=old_t+dt_2;
x_hit_2=(x_2 >=-length/2 && x_2<=length/2);
}

if (prop_case>8)
{
x_3=old_x+vx*dt_3;
y_3=old_y+vy*dt_3;
z_3=old_z+vz*dt_3;
t_3=old_t+dt_3;
x_hit_3=(x_3 >=-length/2 && x_3<=length/2);
}

//printf("x_hit=%d, x_hit_2=%d, x_hit_3=%d\n",x_hit, x_hit_2, x_hit_3);
//printf("dt=%f, dt_2=%f, dt_3=%f\n",dt,dt_2,dt_3);
// printf("x=%f,y=%f,z=%f,vx=%f,vy=%f,vz=%f\n",x,y,z,vx,vy,vz);

// printf("x=%f, length/2=%f\n",x, length/2);


if (x_hit || x_hit_2 || x_hit_3){
//if (x >=-length/2 && x<=length/2){ //check if neutron is within x limits of the mirror. If so, check if it is within y limits.


//define the ellipse
b=smallaxis/2;

f=(focus_e-focus_s)*0.5;

 asquared=f*f+b*b;
 a=sqrt(asquared);

xprime=-f-focus_s+mirror_start+length/2+x; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse

//ymirror=b*sqrt(1-xprime*xprime/(f*f)); //following Kaspars convention, assuming f~=a (valid for most elliptic guides normally used)

ymirror=b*sqrt(1-xprime*xprime/asquared);



xprime_2=-f-focus_s+mirror_start+length/2+x_2; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse
ymirror_2=b*sqrt(1-xprime_2*xprime_2/asquared);

xprime_3=-f-focus_s+mirror_start+length/2+x_3; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse
ymirror_3=b*sqrt(1-xprime_3*xprime_3/asquared);

if (guide_start>mirror_start){ //If (part of the) mirror is outside the guide, the mirror can be extended
if (  x<-length/2+guide_start-mirror_start) {
ymirror=yheight/2;
}

if (  x_2<-length/2+guide_start-mirror_start) {
ymirror_2=yheight/2;
}

if (  x_3<-length/2+guide_start-mirror_start) {
ymirror_3=yheight/2;
}




}











// printf("ymirror=%f, y=%f\n",ymirror, y);
intersect = ( y>=-ymirror && y<=ymirror && x >=-length/2 && x<=length/2);

if (prop_case>0) {
intersect_2 = ( y_2>=-ymirror && y_2<=ymirror && x_2 >=-length/2 && x_2<=length/2);
}
if (prop_case>8){
intersect_3 = ( y_3>=-ymirror && y_3<=ymirror && x_3 >=-length/2 && x_3<=length/2);
}

//printf("y_2=%f, ymirror=%f\n",y_2,ymirror);

//printf("\nintersect=%d, intersect_2=%d, intersect_3=%d, prop_case=%d\n",intersect, intersect_2, intersect_3, prop_case);

//printf("x=%f,y=%f,z=%f,t=%f\n",x,y,z,t);
//printf("x_2=%f,y_2=%f,z_2=%f,t_2=%f\n",x_2,y_2,z_2,t_2);
//printf("x_3=%f,y_3=%f,z_3=%f,t_3=%f\n",x_3,y_3,z_3,t_3);

if (!intersect){
if (!intersect_2){
intersect=intersect_3;
x=x_3;
y=y_3;
z=z_3;
t=t_3;
}else{
intersect=intersect_2;
x=x_2;
y=y_2;
z=z_2;
t=t_2;
}
}

//printf("intersect=%d, intersect_2=%d, intersect_3=%d, prop_case=%d\n\n",intersect, intersect_2, intersect_3, prop_case);
//printf("x=%f,y=%f,z=%f,t=%f\n",x,y,z,t);

//printf("z=%f, zcalc=%f\n",z,a1*x*x*x+b1*x*x+c1*x);
//printf("z=%f, zcalc=%f\n",z_2,a1*x_2*x_2*x_2+b1*x_2*x_2+c1*x_2);
//printf("z=%f, zcalc=%f\n",z_3,a1*x_3*x_3*x_3+b1*x_3*x_3+c1*x_3);

    if (intersect) { //if neutron is within ylimits of the mirror handle reflection/transmission

//first find the angle of the mirror. It is given by theta(x)=alpha*x^2+beta*x+gamma1, with theta(-l/2)=theta1, theta(0)=theta2, theta(l/2)=theta3

alpha1=2*(theta_1+theta_3-2*theta_2)/(length*length);
beta1=(theta_3-theta_1)/length;
gamma1=theta_2;

theta_m=alpha1*x*x+beta1*x+gamma1; // angle of mirror.

//The vector normal to the mirror is e_n= sin(theta)*e_x-cos(theta)*e_z

//find amplitude of v in direction of e_n:

sin_theta_m=sin(theta_m*DEG2RAD);
cos_theta_m=cos(theta_m*DEG2RAD);

v_n=sin_theta_m*vx-cos_theta_m*vz;


q=fabs(2.0*v_n*V2Q);

double R0=0.99;
double Qc=0.0217;
double m_value=m*0.9853+0.1978;
double W=-0.0002*m_value+0.0022;
double alpha=0.1204*m_value+5.0944;
double beta=-7.6251*m_value+68.1137;

if (m_value<=3)
{alpha=m_value;
beta=0;}





      /* Reflectivity (see component Guide). */
      if(m == 0)
        ABSORB;
      if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0"))
         Ref=Table_Value(pTable, q, 1);
      else {
          Ref = R0;
          if(q > Qc)
          {
            double arg = (q-m_value*Qc)/W;
            if(arg < 10)
              Ref *= .5*(1-tanh(arg))*(1-alpha*(q-Qc)+beta*(q-Qc)*(q-Qc)); //matches data from Swiss Neutronics
            else  Ref=0;
          }
      }
      if (Ref < 0) Ref=0;
      else if (Ref > 1) Ref=1;


//Now comes actual reflection/transmission
      if (!transmit) { //all neutrons are reflected
        if (!Ref) ABSORB;
        p *= Ref;

//handle reflection: change v_n -->-v_n

vx=old_vx*(cos_theta_m*cos_theta_m-sin_theta_m*sin_theta_m)+old_vz*(2*cos_theta_m*sin_theta_m);
vz=old_vx*(2*cos_theta_m*sin_theta_m)+old_vz*(sin_theta_m*sin_theta_m-cos_theta_m*cos_theta_m);

// printf("theta_m=%f, sin_theta_m=%f, cos_theta_m=%f, v_n=%f, old_vx=%f, vx=%f, old_vz=%f, vz=%f\n\n", theta_m, sin_theta_m, cos_theta_m, v_n, old_vx, vx, old_vz, vz);


        SCATTER; 
//printf("line 471.In mirror: x=%f,y=%f,z=%f,t=%f\n",x,y,z,t);
//printf("In mirror: old_vx=%f,old_vy=%f,old_vz=%f,vx=%f,vy=%f,vz=%f,v_n=%f\n",old_vx,old_vy,old_vz,vx,vy,vz,v_n);

      } else { //if neutrons can be transmitted



//calculate absorption.
// substrate
double lambda=(2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
double sin_theta=lambda*q/(4*PI);

//double substrate_path_length=substrate_thickness/sin_theta;
//double coating_path_length=coating_thickness/sin_theta;

double sin_theta_c=Qc/(4*PI);

double theta_diff;
double substrate_path_length;
double coating_path_length;

double remaining_length_through_mirror;

int hit_back_mirror;

if (v_n>0) {
hit_back_mirror=1;} else{
hit_back_mirror=0;}

remaining_length_through_mirror=length/2-x;


if (sin_theta>sin_theta_c*lambda) {
theta_diff=sqrt(sin_theta*sin_theta-sin_theta_c*sin_theta_c*lambda*lambda);
coating_path_length=coating_thickness/theta_diff;
substrate_path_length=substrate_thickness/theta_diff;

	if (coating_path_length>remaining_length_through_mirror){
coating_path_length=remaining_length_through_mirror;
substrate_path_length=0; 
}

	if (substrate_path_length>remaining_length_through_mirror){
substrate_path_length=remaining_length_through_mirror; 
}












} else{

if (hit_back_mirror==0){ //neutron comes from front of mirror
substrate_path_length=0;
coating_path_length=remaining_length_through_mirror;
}else {//neutron comes from behind mirror

substrate_path_length=remaining_length_through_mirror;
coating_path_length=0;
}
}


double mu_substrate=0.0318/lambda+0.0055*lambda-0.0050; //unit: cm^-1
mu_substrate=mu_substrate*100; //unit: m^-1;

//For nickel and titanium coating, the following formular is used:
// mu = rho/m(atom)*sigma_a,thermal*lambda/lambda_thermal

// lambda_thermal=1.798 Å

// rho_nickel=8.908g/cm^3
// m(atom)_nickel=58.6934*1.661*10^-27 kg
// sigma_a,thermal_nickel=4.49*10^-28 m^2

// rho_titanium=4.506g/cm^3
// m(atom)_titanium=47.867*1.661*10^-27 kg
// sigma_a,thermal_titanium=6.09*10^-28 m^2

double Ni_coefficient=22.8180; 
double Ti_coefficient=19.1961;

double mu_coating=(0.5*Ni_coefficient+0.5*Ti_coefficient)*lambda; //it is roughly 50% nickel and 50% titanium



        // transmit when rand > R
        if (Ref == 0 || rand01() >= Ref) { //transmit
if (substrate_thickness>0){ p=p*exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length); //reduce weight of neutrons due to attenuation in the mirror
//x+=(coating_path_length+substrate_path_length)-(coating_thickness+substrate_thickness)/sin_theta;
//printf("xshift is %f \n",(coating_path_length+substrate_path_length)-(coating_thickness+substrate_thickness)/sin_theta);
} 
// printf("line 380\n");
/*
if (v_n>0) {
printf("neutron is transmitted from back of mirror. %f\n",exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length));
}else{
printf("neutron is transmitted from front of mirror. %f\n",exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length));
}
*/
} else {//neutron is reflected
		if (v_n>0 && substrate_thickness>0) { //if neutron comes from behind the mirror
// printf("neutron is reflected from back of mirror. %f\n",Ref*exp(-2*mu_substrate*substrate_path_length-2*mu_coating*coating_path_length));
			p=p*exp(-2*mu_substrate*substrate_path_length-2*mu_coating*coating_path_length);} //else{ //reduce weight of neutrons due to attenuation in the mirror
 // printf("neutron is reflected from front of mirror. %f\n", Ref);}
//handle reflection: change v_n -->-v_n
vx=old_vx*(cos_theta_m*cos_theta_m-sin_theta_m*sin_theta_m)+old_vz*(2*cos_theta_m*sin_theta_m);
vz=old_vx*(2*cos_theta_m*sin_theta_m)+old_vz*(sin_theta_m*sin_theta_m-cos_theta_m*cos_theta_m);
// printf("line 388\n");

}

// printf("theta_m=%f, sin_theta_m=%f, cos_theta_m=%f, v_n=%f, old_vx=%f, vx=%f, old_vz=%f, vz=%f\n\n", theta_m, sin_theta_m, cos_theta_m, v_n, old_vx, vx, old_vz, vz);

//printf("vxvx+vzvz=%f, oldvxoldvx+oldvzoldvz=%f", vx*vx+vz*vz, old_vx*old_vx+old_vz*old_vz);

        SCATTER; 
//printf("line 524.In mirror: x=%f,y=%f,z=%f,t=%f\n",x,y,z,t);
//printf("old_vx=%f,old_vy=%f,old_vz=%f,vx=%f,vy=%f,vz=%f,v_n=%f\n",old_vx,old_vy,old_vz,vx,vy,vz,v_n);
//after transmission or reflection
      } //end } else { after if (!transmit) {
    } 


 


} // end if (x >=-length/2 && x<=length/2)

// printf("intersect=%d\n",intersect);

   if (!intersect) {
      /* No intersection: restore neutron position. */
      x = old_x;
      y = old_y;

      z = old_z;
      t = old_t;
// printf("line 409\n");

    }
  

} //end if (dt>0)


%}

MCDISPLAY
%{



/*
if (xcenter==0){

xstart=0;
xend=length;
xprime_start=-a+mirror_start+xstart;
ystart=b*sqrt(1-xprime_start*xprime_start/asquared);

xprime_end=-a+mirror_start+xend;
yend=b*sqrt(1-xprime_end*xprime_end/asquared);

}
*/

/*
if (xcenter==1){
xstart=-length/2;
xend=length/2;


xprime_start=-a+mirror_start+xstart+length/2;
ystart=b*sqrt(1-xprime_start*xprime_start/asquared);

xprime_end=-a+mirror_start+xend+length/2;
yend=b*sqrt(1-xprime_end*xprime_end/asquared);
}

line(xstart,-ystart,0,xstart,ystart,0);
line(xend,-yend,0,xend,yend,0);
line(xstart,-ystart,0,xend,-yend,0);
line(xstart,ystart,0,xend,yend,0);
*/

double xstart;
double xend;
double xprime_start;
double ystart;

double xprime_end;
double yend;

double focus_2;
double focus_1;
double b;
double f;
double asquared;
double a;


int n_lines;
int j=0;
double xprimepos[51];
double ypos[51];
double x_plot[51];
double zpos[51];
double xstep;



focus_2=focus_e-mirror_start; //focus in local coordinates
focus_1=focus_s-mirror_start;

b=smallaxis/2;

f=(focus_2-focus_1)*0.5;
 asquared=f*f+b*b;
 a=sqrt(asquared);


xstart=-length/2;
xend=length/2;

n_lines=50;

xstep=length/n_lines;







double a1, b1, c1;
double tan_theta_1;
double tan_theta_2;
double tan_theta_3;



//mirror is defined by z(x)=a1x^3+b1x^2+c1x+d1, with dz/dx|x=-length/2=tan(theta_1), dz/dx|x=0=tan(theta_2), dz/dx|x=length/2=tan(theta3), z(0)=0. (d1=0)

tan_theta_1=tan(theta_1*DEG2RAD);
tan_theta_2=tan(theta_2*DEG2RAD);
tan_theta_3=tan(theta_3*DEG2RAD);


a1=2.0/3.0*(tan_theta_1+tan_theta_3-2.0*tan_theta_2)/(length*length);
b1=(tan_theta_3-tan_theta_1)/(2.0*length);
c1=tan_theta_2;



for (j=0; j<n_lines+1; j++)
{
xprimepos[j]=-f-focus_s+mirror_start+length/2+xstart+xstep*j;

ypos[j]=b*sqrt(1-xprimepos[j]*xprimepos[j]/asquared); //correct

if (guide_start>mirror_start){
if (  xstart+xstep*j<-length/2+guide_start-mirror_start) {
ypos[j]=yheight/2;
}
}



// ypos[j]=b*sqrt(1-xprimepos[j]*xprimepos[j]/(f*f)); //following convention in Kaspar's elliptic guide..
// printf("xprimepos[j]=%f,f*f=%f, ypos[j]=%f\n",xprimepos[j],f*f,ypos[j]);

x_plot[j]=xstart+xstep*j;
zpos[j]=a1*x_plot[j]*x_plot[j]*x_plot[j]+b1*x_plot[j]*x_plot[j]+c1*x_plot[j];
}

for (j=0; j<n_lines; j++)
{
line(x_plot[j], -ypos[j], zpos[j], x_plot[j+1], -ypos[j+1],zpos[j+1]);
line(x_plot[j], ypos[j], zpos[j], x_plot[j+1], ypos[j+1],zpos[j+1]);
}


line(x_plot[0],-ypos[0],zpos[0],x_plot[0],ypos[0],zpos[0]);
line(x_plot[50],-ypos[50],zpos[50],x_plot[50],ypos[50],zpos[50]);




//printf("ypos0=%f xpos0=%f ypos50=%f xpos50=%f",ypos[0], x_plot[0], ypos[50], x_plot[50]);

/*  double xmax, xmin, ymax, ymin;
  

  if (center == 0) {
    xmax= x1; xmin=0;
    ymax= yheight; ymin=0;
  } else {
    xmax= x1/2; xmin=-xmax;
    ymax= yheight/2; ymin=-ymax;
  }
  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