ESS-Estia/Estia-McStas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
c85792ef22c77cd7d0cacee01d0099cfbfc8f073
Estia-McStas/simulation/Polariser.comp
T

1263 lines
56 KiB
Plaintext
Raw Blame History

/**********************************************************************************************************************
* McStas instrument definitions: http://www.mcstas.org
*
* Instrument: ESS Reflectometer Estia (ESS_Reflectometer_Estia)
* Component: Polariser
* Authors: Developed at the Laboratory for Neutron Scattering and Imaging (LNS) by Evangelos Matsinos and Artur Glavic;
* for information and questions, e-mail artur.glavic@psi.ch
* Origin: Paul Scherrer Institut (PSI), Villigen, Switzerland
* Release: McStas x.x
* Version: 1.0
* Date: 17/03/2018
*
* Description
* -----------
* The module provides a simulation of the neutron-beam propagation through the polarisers. Polarisers are selectors of
* one spin state (the neutrons with opposite spin are mostly reflected away from the beamline). Polarisers are usually
* silicon films of a specific geometry - spiral (y,z) cross section - and are appropriately coated in order to enable
* the transmission of the incident neutrons of one spin state.
*
* Input parameters
* ----------------
* enable_ref: [] If 0, the material of the polariser does not refract (to enable testing)
* abs_out: [] If 0, the neutrons, which do not intersect the mirror, are propagated (i.e., not absorbed)
* abs_single: [] If 0, the neutrons, which do only intersect the first mirror, are propagated (i.e., not absorbed)
* abs_ref: [] If 0, the neutrons, which are reflected from the mirror, are propagated (i.e., not absorbed)
* both_coated: [] If 0, only the up-stream surface is coated
* lin: [m] Entrance point (own CS): if negative, proximal polariser; if positive, distal polariser
* length: [m] Length (positive definite) of the (proximal or distal) polariser in the z direction (own CS)
* delta_theta: [rad] The angle subtended by the polariser at the focal point
* h1: [m] The total length of the polariser in the x direction at z=dZmin (own CS)
* h2: [m] The total length of the polariser in the x direction at z=dZmax (own CS)
* R0: [] The low-angle reflectivity (StdReflecFunc)
* alpha: [A] The slope of the reflectivity for the spin-up state (StdReflecFunc)
* m_d: [] The m-value of the coating for the spin-down state (StdReflecFunc)
* m_u: [] The m-value of the coating for the spin-up state (StdReflecFunc)
* m_substrate: [] The m-value of the material of the polariser
* d_substrate: [m] The thickness of the polariser (typically, a fraction of 1mm)
* m_residual: [] The m-value used for residual Q^-4 dependent reflectivity
* W: [] The width of the supermirror cut-off (StdReflecFunc)
* T_loss: [] Amount of intensity lost due to transmission.
* reflect_d: [str] The file path to the reflectivity R for the spin-down state. Format: q [A^-1] R
* reflect_u: [str] The file path to the reflectivity R for the spin-up state. Format: q [A^-1] R
*
* Output
* ------
*
*
**********************************************************************************************************************/
DEFINE COMPONENT Polariser
DEFINITION PARAMETERS ()
SETTING PARAMETERS (int enable_ref = 1,
int abs_out = 1,
int abs_single = 0,
int abs_ref = 1,
int both_coated = 1,
lin = 170.0e-03,
length = 225.0e-03,
delta_theta = 0.0314159,
h1 = 300.0e-03,
h2 = 400.0e-03,
R0 = 1.0,
alpha = 6.49,
m_d = 0.5,
m_u = 4.0,
m_substrate = 0.469522,
d_substrate = 0.5e-03,
m_residual = 0.57,
W = 0.0033333,
T_loss=4.0e3,
string reflect_d = "ReflectivitySpinDown.dat",
string reflect_u = "ReflectivitySpinUp.dat")
/*
OUTPUT PARAMETERS (dDThetaHalf,
dZmin,
dZmax,
ab,
bb,
dOffsetY,
qc_substrate,
nRflctUp,
nRflctDn,
dGAx,
dGAy,
dGAz)
DECLARE
%{
%}
*/
SHARE
%{
%include "pol-lib"
%include "read_table-lib"
%include "ref-lib"
/**********************************************************************************************************************
Declaration of the variables used in the McStas method Polariser */
t_Table* pTableDn=NULL; /* The table for the spin-down reflectivity */
t_Table* pTableUp=NULL; /* The table for the spin-up reflectivity */
int debug_print=0; /* For debugging purposes */
int abs_reason[]={0,0,0,0,0,0,0,0}; /* Auxiliary memory used for flagging the absorbed neutrons */
/* abs_reason[0]: unable to obtain intersection with the first surface of the polariser */
/* abs_reason[1]: unreliable intersection with the first surface of the polariser */
/* abs_reason[2]: out of bounds (in the x direction) intersection with the first surface of the polariser */
/* abs_reason[3]: reflection at the intersection with the first surface of the polariser */
/* abs_reason[4]: unable to obtain intersection with the second surface of the polariser */
/* abs_reason[5]: unreliable intersection with the second surface of the polariser */
/* abs_reason[6]: out of bounds (in the x direction) intersection with the second surface of the polariser */
/* abs_reason[7]: unphysical (null velocity, erroneous refractive index, negative propagation time, non-positive weight, etc.) */
double pdTmp[]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}; /* Auxiliary memory */
int nRflctDn,nRflctUp;
double dTlrnce,dDThetaHalf,qc_Ni,qc_substrate,ab,bb,dOffsetY,dZmin,dZmax,dGAx,dGAy,dGAz;
double dSpdLght,dMassNeut,dMassNeut,dhbarc;
Coords localG; /* Gravitational-acceleration vector */
Coords eCRD_3DU1; /* Constant unit vector (x direction) */
/*
-----------------------------------------------------------------------------------------------------------------------
Definition of functions used in the main body of the method
*/
double FTN1(double dX) {double dTmp=exp(pdTmp[10]*dX); double dSin=sin(dX);
return dTmp*((pdTmp[0]*dTmp*dSin+pdTmp[1])*dSin+pdTmp[2]*cos(dX))+pdTmp[3];};
double FTN2(double dX) {double dTmp=exp(pdTmp[10]*dX); double dCos=cos(dX);
return dTmp*((pdTmp[0]*dTmp*dCos+pdTmp[1])*dCos+pdTmp[2]*sin(dX))+pdTmp[3];};
double FTN3(double dX) {double dT=(exp(pdTmp[10]*dX)*(pdTmp[0]*sin(dX)+pdTmp[1]*cos(dX))+pdTmp[2])/pdTmp[3];
return pdTmp[4]+(pdTmp[5]+pdTmp[6]*dT)*dT+(pdTmp[7]+(pdTmp[8]+pdTmp[9]*dT)*dT)*tan(dX);};
double FindRoot(double (*FUNCTION)(double ),const double dX1,const double dX2,const double dAccuracy,int* pnSuccess);
/*
-----------------------------------------------------------------------------------------------------------------------
Given a function FUNCTION defined in the interval [dX1,dX2], subdivide the interval into nNoInts equally spaced
segments, and search for zero crossings of the function. The value of pnNoRoots at input is the maximal number of roots
sought, and is at output reset to the number of bracketing pairs pdX1[1,...,pnNoRoots], pdX2[1,...,pnNoRoots] found.
*/
void BracketRoots(double (*FUNCTION)(double ),
const double dX1,
const double dX2,
const int nNoInts,
double* pdX1,
double* pdX2,
int* pnNoRoots) {
int nbb,nCnt;
double dX,dValue1,dValue2,dStep;
nbb=0;
dStep=(dX2-dX1)/(double)nNoInts; /* Determine the spacing appropriate to the mesh */
dValue1=(*FUNCTION)(dX=dX1);
for(nCnt=0;nCnt<nNoInts;++nCnt) { /* Loop over all intervals */
dValue2=(*FUNCTION)(dX += dStep);
if(dValue2*dValue1 <= 0.0) {
pdX1[nbb]=dX-dStep;
pdX2[nbb++]=dX;
if(*pnNoRoots == nbb)return;
}
dValue1=dValue2;
}
*pnNoRoots=nbb;
}
/*
-----------------------------------------------------------------------------------------------------------------------
Given a function FUNCTION and an initial guess range from dX1 to dX2, the method expands the range geometrically until
a root is bracketed by the returned values dX1 and dX2 (in which case the method returns 1) or until the range becomes
unacceptably large (in which case the method returns 0).
*/
int BracketRoot(double (*FUNCTION)(double ),
double* pdX1,
double* pdX2) {
const int nNoIters=100;
const double dFactor=1.6;
int nCnt;
double dValue1,dValue2;
if(*pdX1 == *pdX2)*pdX2 += (double)nNoIters;
dValue1=(*FUNCTION)(*pdX1);
dValue2=(*FUNCTION)(*pdX2);
for(nCnt=0;nCnt<nNoIters;++nCnt) {
if(dValue1*dValue2<0.0)return 1;
if(fabs(dValue1)<fabs(dValue2)) {
dValue1=(*FUNCTION)(*pdX1 += dFactor*(*pdX1-*pdX2));
} else {
dValue2=(*FUNCTION)(*pdX2 += dFactor*(*pdX2-*pdX1));
}
}
return 0;
}
/*
-----------------------------------------------------------------------------------------------------------------------
Assertion of the appropriate intersection with the spiral surface of the distal polariser.
*/
int IntersectsDistalPolariser(const double dX,
const double dY,
const double dZ,
const double dH1In,
const double dH2In,
const double dZminIn,
const double dZmaxIn) {
const double dHlfSum=0.5*(dZminIn+dZmaxIn);
const double dHlfDff=0.5*(dZmaxIn-dZminIn);
if(fabs(dZ-dHlfSum)>dHlfDff)return 0;
if(fabs(dX)>0.5*((dH2In-dH1In)*fabs(dZ)+dZmaxIn*dH1In-dZminIn*dH2In)/(dZmaxIn-dZminIn))return 0;
return 1;
}
/*
-----------------------------------------------------------------------------------------------------------------------
Assertion of the appropriate intersection with the spiral surface of the proximal polariser.
*/
int IntersectsProximalPolariser(const double dX,
const double dY,
const double dZ,
const double dH1In,
const double dH2In,
const double dZminIn,
const double dZmaxIn) {
const double dHlfSum=0.5*(dZminIn+dZmaxIn);
const double dHlfDff=0.5*(dZmaxIn-dZminIn);
if(fabs(dZ+dHlfSum)>dHlfDff)return 0;
if(fabs(dX)>0.5*((dH2In-dH1In)*fabs(dZ)+dZmaxIn*dH1In-dZminIn*dH2In)/(dZmaxIn-dZminIn))return 0;
return 1;
}
/*
-----------------------------------------------------------------------------------------------------------------------
The method checks whether the input vector is null
*/
int IsVectorNull(const Coords eCRD_3D) {
const double dTlrnce_IsVectorNull=1.0e-08; /* The default tolerance level/inaccuracy in the estimates */
double dX,dY,dZ;
coords_get(eCRD_3D,&dX,&dY,&dZ);
if((fabs(dX)>dTlrnce_IsVectorNull) || (fabs(dY)>dTlrnce_IsVectorNull) || (fabs(dZ)>dTlrnce_IsVectorNull))return 0;
return 1;
}
/*
-----------------------------------------------------------------------------------------------------------------------
The method returns the magnitude of a given input vector
*/
int ObtainMagnitudeOfVector(const Coords eCRD_3D,
double* dLength) {
double dX,dY,dZ;
*dLength=0.0;
coords_get(eCRD_3D,&dX,&dY,&dZ);
*dLength=dX*dX+dY*dY+dZ*dZ;
*dLength=(*dLength>0.0) ? sqrt(*dLength) : 0.0;
return (*dLength>0.0) ? 1 : 0;
}
/*
-----------------------------------------------------------------------------------------------------------------------
Using the van Wijngaarden-Dekker-Brent method, the module returns the root of a function FUNCTION. A root is assumed to
exist in the [dX1,dX2] interval. The returned value is refined until the accuracy of the estimate is at least equal to
dAccuracy.
*/
#define SHFT_FindRoot(dA,dB,dC) dA=dB;dB=dC;dC=dA;
#define SIGN_FindRoot(dA,dB) ((dB) >= 0.0 ? fabs(dA) : -fabs(dA))
double FindRoot(double (*FUNCTION)(double ),
const double dX1,
const double dX2,
const double dAccuracy,
int* pnSuccess) {
const int nNoIters=100;
const double dEps=1.0e-10;
int nIter,nNoRoots;
double dValue3,dMin1,dMin2,dP,dQ,dR,dS,dTmp,dX;
double pdX1[nNoIters],pdX2[nNoIters];
double dA=dX1;
double dB=dX2;
double dC=dX2;
double dD=0.0;
double dE=0.0;
double dValue1=(*FUNCTION)(dA);
double dValue2=(*FUNCTION)(dB);
*pnSuccess=0;
if(dValue1*dValue2>0.0) { /* Attempt to bracket the minimum */
if(BracketRoot(FUNCTION,&dA,&dB)) {
dC=dB;
dValue1=(*FUNCTION)(dA);
dValue2=(*FUNCTION)(dB);
} else {
BracketRoots(FUNCTION,dX1,dX2,nNoIters,pdX1,pdX2,&nNoRoots);
if(nNoRoots) { /* First root to be considered */
dA=pdX1[0];
dB=dC=pdX2[0];
dValue1=(*FUNCTION)(dA);
dValue2=(*FUNCTION)(dB);
} else {
return 0.0;
}
}
}
dValue3=dValue2;
for(nIter=0;nIter<nNoIters;++nIter) {
if((dValue2>0.0 && dValue3>0.0) || (dValue2<0.0 && dValue3<0.0)) {
dC=dA;
dValue3=dValue1;
dE=dD=dB-dA;
}
if(fabs(dValue3)<fabs(dValue2)) {
SHFT_FindRoot(dA,dB,dC)
SHFT_FindRoot(dValue1,dValue2,dValue3)
}
dTmp=2.0*dEps*fabs(dB)+0.5*dAccuracy;
dX=0.5*(dC-dB);
if((fabs(dX) <= dTmp) || (dValue2 == 0.0)) {
*pnSuccess=1;
return dB;
}
if((fabs(dE) >= dTmp) && (fabs(dValue1)>fabs(dValue2))) {
dS=dValue2/dValue1;
if(dA == dC) {
dP=2.0*dX*dS;
dQ=1.0-dS;
} else {
dQ=dValue1/dValue3;
dR=dValue2/dValue3;
dP=dS*(2.0*dX*dQ*(dQ-dR)-(dB-dA)*(dR-1.0));
dQ=(dQ-1.0)*(dR-1.0)*(dS-1.0);
}
if(dP>0.0)dQ=-dQ;
dP=fabs(dP);
dMin1=3.0*dX*dQ-fabs(dTmp*dQ);
dMin2=fabs(dE*dQ);
if(2.0*dP<(dMin1<dMin2 ? dMin1 : dMin2)) {
dE=dD;
dD=dP/dQ;
} else {
dD=dX;
dE=dD;
}
} else {
dD=dX;
dE=dD;
}
dA=dB;
dValue1=dValue2;
if(fabs(dD)>dTmp) {
dB += dD;
} else {
dB += SIGN_FindRoot(dTmp,dX);
}
dValue2=(*FUNCTION)(dB);
}
return 0.0;
}
#undef SHFT_FindRoot
#undef SIGN_FindRoot
/*
-----------------------------------------------------------------------------------------------------------------------
The method returns the refractive index
*/
double ObtainRefractiveIndex(const double dWL, /* Wavelength (A) */
const double dQc) { /* The critical q value (A^-1) */
double dRI=0.25*dWL*dQc/PI;
return (dRI<1.0) ? sqrt(1.0-dRI*dRI) : 0.0;
}
/*
-----------------------------------------------------------------------------------------------------------------------
The method returns the spin-dependent reflection weight
*/
double ObtainReflectivity(const int nRflct, /* If 0, StdReflecFunc is used; otherwise, the supplied table is used */
const double dVT, /* Velocity perpendicular to the surface [m s^-1] */
const double dR0, /* Used in StdReflecFunc */
const double dQc, /* Used in StdReflecFunc */
const double dSlp, /* Used in StdReflecFunc */
const double dMctg, /* Used in StdReflecFunc */
const double dCutOff, /* Used in StdReflecFunc */
const double m_residual,
t_Table* pt_Table) { /* Used in TableReflecFunc */
double dW;
double minW;
double dqT=fabs(2.0*dVT*V2Q); /* V2Q is equal to the ratio of the neutron mass to Planck's reduced constant ~ 1.58825e-03 [A^-1 (m s^-1)^-1] */
if(!nRflct) {
double pPar[]={dR0,dQc,dSlp,dMctg,dCutOff};
StdReflecFunc(dqT,pPar,&dW);
if ((dqT/dQc)>dMctg) {
minW=pow((dqT/dQc),-4)/pow(m_residual, -4);
dW=fmax(minW,dW);
}
} else {
TableReflecFunc(dqT,pt_Table,&dW);
}
if(dW<0.0)dW=0.0;
if(dW>1.0)dW=1.0;
return dW;
}
/*
-----------------------------------------------------------------------------------------------------------------------
The method determines the characteristics of the intersection of the neutron with the distal part of the polariser.
When changing the McStas variables, characterising each neutron, without using the dedicated McStas functions, a
suspicious behaviour of the application was found; this shortcoming is probably due to the direct access of the McStas
allocated memory. To remedy this, the position and velocity components of each neutron are passed as input (pdNCt) to
this module, and the neutron state is updated by using only the PROP_DT function.
*/
void DetermineIntersectionWithDistalPolariser(double* pdNCt, /* INPUT: The position and velocity vectors of the neutron */
double* pdGA, /* INPUT: Gravitational acceleration in the y and z directions */
double* pdC, /* INPUT: Various constants */
/* OUTPUT FROM HERE ON */
int* pnType, /* 0: regular neutrons: well-defined intersections */
/* 1: irregular neutrons to be absorbed (or propagated, if the method is called for the second surface and abs_single is set to 0) (failed cases) */
/* 2: irregular neutrons to be propagated (failed cases) */
double* pdRt, /* The angle theta corresponding to the intersection (only for regular neutrons) */
double* pdT) { /* The time needed for the neutron to reach the intersection (only for regular neutrons) */
int nSuccess;
/* Starting position and velocity vectors of the neutron */
double dX=pdNCt[0]; /* Neutron x coordinate (own CS) */
double dY=pdNCt[1]; /* Neutron y coordinate (own CS) */
double dZ=pdNCt[2]; /* Neutron z coordinate (own CS) */
double dVx=pdNCt[3]; /* Neutron x velocity (own CS) */
double dVy=pdNCt[4]; /* Neutron y velocity (own CS) */
double dVz=pdNCt[5]; /* Neutron z velocity (own CS) */
/* Gravitational-acceleration components */
double dGAy=pdGA[0]; /* Gravitational acceleration in the y direction */
double dGAz=pdGA[1]; /* Gravitational acceleration in the z direction */
/* Various constants */
double dOffY=pdC[0]; /* Offset in y, to account for the two surfaces of the polariser */
double dHalfRange=pdC[1]; /* The limit of the angle theta */
double dTol=pdC[2]; /* The tolerance used in the estimation of the roots */
int nAbsIrr=(int)(pdC[3]); /* Input to the Polariser component */
int nOff=(fabs(dOffY)>dTol) ? 0 : 4;
/* Initialisation */
*pnType=0; /* All neutrons are assumed regular in the beginning */
*pdRt=*pdT=0.0;
/* Determination of the intersection (corresponding angle theta) and of the time required for the neutron to reach it */
if(dGAy == 0.0) { /* Case: no gravity in the y direction (simplified case) */
pdTmp[0]=0.5*dGAz*ab*ab;
pdTmp[1]=ab*(dVy*dVz-dGAz*(dY-dOffY));
pdTmp[2]=-ab*dVy*dVy;
pdTmp[3]=dVy*dVy*dZ-dVy*dVz*(dY-dOffY)+0.5*dGAz*(dY-dOffY)*(dY-dOffY);
if(FTN1(-dHalfRange)*FTN1(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN1,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdT=(ab*exp(bb*(*pdRt))*sin(*pdRt)-(dY-dOffY))/dVy; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
} else if(dGAz == 0.0) { /* Case: no gravity in the z direction (simplified case) */
pdTmp[0]=0.5*dGAy*ab*ab;
pdTmp[1]=ab*(dVy*dVz-dGAy*dZ);
pdTmp[2]=-ab*dVz*dVz;
pdTmp[3]=dVz*dVz*(dY-dOffY)-dVy*dVz*dZ+0.5*dGAy*dZ*dZ;
if(FTN2(-dHalfRange)*FTN2(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN2,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdT=(ab*exp(bb*(*pdRt))*cos(*pdRt)-dZ)/dVz; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
} else { /* Case: gravitational acceleration components in the y and z directions (general case) */
pdTmp[0]=ab*dGAz;
pdTmp[1]=-ab*dGAy;
pdTmp[2]=dGAy*dZ-dGAz*(dY-dOffY);
pdTmp[3]=dGAz*dVy-dGAy*dVz;
pdTmp[4]=dY-dOffY;
pdTmp[5]=dVy;
pdTmp[6]=0.5*dGAy;
pdTmp[7]=-dZ;
pdTmp[8]=-dVz;
pdTmp[9]=-0.5*dGAz;
if(FTN3(-dHalfRange)*FTN3(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN3,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdT=(exp(bb*(*pdRt))*(pdTmp[0]*sin(*pdRt)+pdTmp[1]*cos(*pdRt))+pdTmp[2])/pdTmp[3]; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
} /* Endif(dGAy == 0.0) */
return;
}
/*
-----------------------------------------------------------------------------------------------------------------------
The method determines the characteristics of the intersection of the neutron with the proximal part of the polariser.
When changing the McStas variables, characterising each neutron, without using the dedicated McStas functions, a
suspicious behaviour of the application was found; this shortcoming is probably due to the direct access of the McStas
allocated memory. To remedy this, the position and velocity components of each neutron are passed as input (pdNCt) to
this module, and the neutron state is updated by using only the PROP_DT function.
*/
void DetermineIntersectionWithProximalPolariser(double* pdNCt, /* INPUT: The position and velocity vectors of the neutron */
double* pdGA, /* INPUT: Gravitational acceleration in the y and z directions */
double* pdC, /* INPUT: Various constants */
/* OUTPUT FROM HERE ON */
int* pnType, /* 0: regular neutrons: well-defined intersections */
/* 1: irregular neutrons to be absorbed (or propagated, if the method is called for the second surface and abs_single is set to 0) (failed cases) */
/* 2: irregular neutrons to be propagated (failed cases) */
double* pdRt, /* The angle theta corresponding to the intersection (only for regular neutrons) */
double* pdT) { /* The time needed for the neutron to reach the intersection (only for regular neutrons) */
int nSuccess;
/* Starting position and velocity vectors of the neutron */
double dX=pdNCt[0]; /* Neutron x coordinate (own CS) */
double dY=pdNCt[1]; /* Neutron y coordinate (own CS) */
double dZ=pdNCt[2]; /* Neutron z coordinate (own CS) */
double dVx=pdNCt[3]; /* Neutron x velocity (own CS) */
double dVy=pdNCt[4]; /* Neutron y velocity (own CS) */
double dVz=pdNCt[5]; /* Neutron z velocity (own CS) */
/* Gravitational-acceleration components */
double dGAy=pdGA[0]; /* Gravitational acceleration in the y direction */
double dGAz=pdGA[1]; /* Gravitational acceleration in the z direction */
/* Various constants */
double dOffY=pdC[0]; /* Offset in y, to account for the two surfaces of the polariser */
double dHalfRange=pdC[1]; /* The limit of the angle theta */
double dTol=pdC[2]; /* The tolerance used in the estimation of the roots */
int nAbsIrr=(int)(pdC[3]); /* Input to the Polariser component */
int nOff=(fabs(dOffY)<dTol) ? 0 : 4;
/* Initialisation */
*pnType=0; /* All neutrons are assumed regular in the beginning */
*pdRt=*pdT=0.0;
/* Determination of the intersection (corresponding angle theta) and of the time required for the neutron to reach it */
if(dGAy == 0.0) { /* Case: no gravity in the y direction (simplified case) */
pdTmp[0]=0.5*dGAz*ab*ab;
pdTmp[1]=ab*(dVy*dVz-dGAz*(dY-dOffY));
pdTmp[2]=ab*dVy*dVy;
pdTmp[3]=dVy*dVy*dZ-dVy*dVz*(dY-dOffY)+0.5*dGAz*(dY-dOffY)*(dY-dOffY);
if(FTN1(-dHalfRange)*FTN1(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN1,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdT=(ab*exp(bb*(*pdRt))*sin(*pdRt)-(dY-dOffY))/dVy; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
} else if(dGAz == 0.0) { /* Case: no gravity in the z direction (simplified case) */
pdTmp[0]=0.5*dGAy*ab*ab;
pdTmp[1]=ab*(dGAy*dZ-dVy*dVz);
pdTmp[2]=-ab*dVz*dVz;
pdTmp[3]=dVz*dVz*(dY-dOffY)-dVy*dVz*dZ+0.5*dGAy*dZ*dZ;
if(FTN2(-dHalfRange)*FTN2(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN2,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdT=-(ab*exp(bb*(*pdRt))*cos(*pdRt)+dZ)/dVz; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
} else { /* Case: gravitational acceleration components in the y and z directions (general case) */
pdTmp[0]=ab*dGAz;
pdTmp[1]=ab*dGAy;
pdTmp[2]=dGAy*dZ-dGAz*(dY-dOffY);
pdTmp[3]=dGAz*dVy-dGAy*dVz;
pdTmp[4]=dY-dOffY;
pdTmp[5]=dVy;
pdTmp[6]=0.5*dGAy;
pdTmp[7]=dZ;
pdTmp[8]=dVz;
pdTmp[9]=0.5*dGAz;
if(FTN3(-dHalfRange)*FTN3(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN3,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdT=(exp(bb*(*pdRt))*(pdTmp[0]*sin(*pdRt)+pdTmp[1]*cos(*pdRt))+pdTmp[2])/pdTmp[3]; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
} /* Endif(dGAy == 0.0) */
return;
}
/*
-----------------------------------------------------------------------------------------------------------------------
The method estimates the time needed for the neutron to reach a specific z
*/
void EstimateExitTime(double* pdT) {
double dTmp;
*pdT=0.0;
if(pdTmp[2] == 0.0) {
if(pdTmp[1] == 0.0)return;
*pdT=-pdTmp[0]/pdTmp[1]; /* Solution under constant velocity */
} else {
dTmp=pdTmp[1]*pdTmp[1]-4.0*pdTmp[0]*pdTmp[2];
if(dTmp<0.0) {
if(pdTmp[1] == 0.0)return;
*pdT=-pdTmp[0]/pdTmp[1]; /* Estimate on the basis of the quadratic form failed: solution under constant velocity */
} else {
dTmp=sqrt(dTmp);
*pdT=0.5*(dTmp-pdTmp[1])/pdTmp[2];
if(*pdT<0.0)*pdT=-pdTmp[0]/pdTmp[1]; /* Estimate on the basis of the quadratic form failed: solution under constant velocity */
} /* Endif(dTmp<0.0) */
} /* Endif(pdTmp[2] == 0.0) */
return;
}
%}
/**********************************************************************************************************************
Initialisation of variables and constructs used in the McStas method Polariser */
INITIALIZE
%{
/* Vefirication of the input */
if((length <= 0.0) || (h1<0.0) || (h2<0.0) || (delta_theta <= 0.0) ||
(d_substrate <= 0.0) || (m_d <= 0.0) || (m_u <= 0.0) || (fabs(R0-0.5)>0.5+dTlrnce))
exit(fprintf(stderr,"%s: Invalid input parameter\n",NAME_CURRENT_COMP));
/* Project-related constants */
dTlrnce=1.0e-06; /* The default tolerance level/inaccuracy in the estimates */
/* Physical constants */
dSpdLght=299792458.0; /* [m s^-1] The speed of light */
dMassNeut=939.5654133; /* [MeV] The neutron rest energy */
dhbarc=197.3269788; /* [MeV fm] The hbar c value */
qc_Ni=sqrt(16.0*PI*9.408e-06); /* [A^-1] The critical q value for Ni; the coherent SLD is used */
/* Derivable constants */
dDThetaHalf=0.5*delta_theta; /* Half of the angle delta_theta; to avoid estimating the half angle at several places */
dZmin=(lin<0.0) ? -(lin+length) : lin;
dZmax=(lin<0.0) ? -lin : lin+length;
ab=sqrt(dZmax*dZmin)/cos(dDThetaHalf); /* [m] The parameter a of the spiral curve */
bb=log(dZmax/dZmin)/delta_theta; /* [rad^-1] The parameter b of the spiral curve */
dOffsetY=ab*sqrt(bb*bb+1.0)*(exp(bb*dDThetaHalf)-exp(-bb*dDThetaHalf))*d_substrate/(bb*(dZmax-dZmin)); /* [m] The separation between the two surfaces of the polariser in the y direction */
qc_substrate=m_substrate*qc_Ni; /* [A^-1] The critical q value for the substrate */
nRflctDn=0; /* If an existing file is given as input for the reflectivity data, this variable will be turned to 1 */
nRflctUp=0; /* If an existing file is given as input for the reflectivity data, this variable will be turned to 1 */
localG=(mcgravitation) ? rot_apply(ROT_A_CURRENT_COMP,coords_set(0.0,-GRAVITY,0.0)) : coords_set(0.0,0.0,0.0);
coords_get(localG,&dGAx,&dGAy,&dGAz);
if(fabs(dGAx)<dTlrnce)dGAx=0.0;
if(fabs(dGAy)<dTlrnce)dGAy=0.0;
if(fabs(dGAz)<dTlrnce)dGAz=0.0;
eCRD_3DU1=coords_set(1.0,0.0,0.0); /* Unit vector in the x direction */
pTableDn=(t_Table*)malloc(sizeof(t_Table));
pTableUp=(t_Table*)malloc(sizeof(t_Table));
if(reflect_d && strlen(reflect_d) && strcmp(reflect_d,"NULL") && strcmp(reflect_d,"0")) {
if(Table_Read(pTableDn,reflect_d,1) <= 0) /* Read the first block data from the input file into pTableDn */
exit(fprintf(stderr,"%s: unable to input the spin-down reflectivity: %s\n",NAME_CURRENT_COMP,reflect_d));
nRflctDn=1;
} /* Endif(reflect_d && strlen(reflect_d) && strcmp(reflect_d,"NULL") && strcmp(reflect_d,"0")) */
if(reflect_u && strlen(reflect_u) && strcmp(reflect_u,"NULL") && strcmp(reflect_u,"0")) {
if(Table_Read(pTableUp,reflect_u,1) <= 0) /* Read the first block data from the input file into pTableUp */
exit(fprintf(stderr,"%s: unable to input the spin-up reflectivity: %s\n",NAME_CURRENT_COMP,reflect_u));
nRflctUp=1;
} /* Endif(reflect_u && strlen(reflect_u) && strcmp(reflect_u,"NULL") && strcmp(reflect_u,"0")) */
if(debug_print) {
fprintf(stdout,"\nInput parameters to %s\n",NAME_CURRENT_COMP);
fprintf(stdout,"==============================\n");
fprintf(stdout,"Inclusion of refraction [] : %d\n",enable_ref);
fprintf(stdout,"Absorption outside polariser [] : %d\n",abs_out);
fprintf(stdout,"Absorption signle interface [] : %d\n",abs_single);
fprintf(stdout,"Absorption of reflected neutrons [] : %d\n",abs_ref);
fprintf(stdout,"Absorption single interface [] : %d\n",abs_single);
fprintf(stdout,"Both sides coated [] : %d\n",both_coated);
fprintf(stdout,"Distance zmin [m] : %f\n",dZmin);
fprintf(stdout,"Distance zmax [m] : %f\n",dZmax);
fprintf(stdout,"Length h1 [m] : %f\n",h1);
fprintf(stdout,"Length h2 [m] : %f\n",h2);
fprintf(stdout,"Angle Dtheta [rad] : %f\n",delta_theta);
fprintf(stdout,"Polariser thickness [m] : %f\n",d_substrate);
fprintf(stdout,"Polariser inner-surface y-offset [m] : %f\n",dOffsetY);
fprintf(stdout,"Low-angle reflectivity [] : %f\n",R0);
fprintf(stdout,"Slope of reflectivity [A^-1] : %f\n",alpha);
fprintf(stdout,"m-value for spin-up [] : %f\n",m_u);
fprintf(stdout,"m-value for spin-down [] : %f\n",m_d);
fprintf(stdout,"m-value of the polariser material [] : %f\n",m_substrate);
fprintf(stdout,"Cut-off width [] : %f\n",W);
fprintf(stdout,"Gravitational acceleration x [m s^-2]: %f\n",dGAx);
fprintf(stdout,"Gravitational acceleration y [m s^-2]: %f\n",dGAy);
fprintf(stdout,"Gravitational acceleration z [m s^-2]: %f\n\n",dGAz);
if(nRflctDn) {
fprintf(stdout,"Spin-down reflectivity file [] : %s\n",reflect_d);
} else {
fprintf(stdout,"Spin-down reflectivity via the function StdReflecFunc\n");
} /* Endif(nRflctDn) */
if(nRflctUp) {
fprintf(stdout,"Spin-up reflectivity file [] : %s\n",reflect_u);
} else {
fprintf(stdout,"Spin-up reflectivity via the function StdReflecFunc\n");
} /* if(nRflctUp) */
if(mcgravitation) {
fprintf(stdout,"Gravitation included\n");
} else {
fprintf(stdout,"Gravitation omitted\n");
} /* Endif(mcgravitation) */
} /* Endif(debug_print) */
%}
/**********************************************************************************************************************
The main body of the method */
TRACE
%{
/* Reinitialize parameters */
dDThetaHalf=0.5*delta_theta; /* Half of the angle delta_theta; to avoid estimating the half angle at several places */
dZmin=(lin<0.0) ? -(lin+length) : lin;
dZmax=(lin<0.0) ? -lin : lin+length;
ab=sqrt(dZmax*dZmin)/cos(dDThetaHalf); /* [m] The parameter a of the spiral curve */
bb=log(dZmax/dZmin)/delta_theta; /* [rad^-1] The parameter b of the spiral curve */
dOffsetY=ab*sqrt(bb*bb+1.0)*(exp(bb*dDThetaHalf)-exp(-bb*dDThetaHalf))*d_substrate/(bb*(dZmax-dZmin)); /* [m] The separation between the two surfaces of the polariser in the y direction */
qc_substrate=m_substrate*qc_Ni; /* [A^-1] The critical q value for the substrate */
int nType,nRefracted,isPolarising;
double dRt,dT,dSlp,dTmp,dTmp1,dTmp2,dTmp3,dX1,dY1,dZ1,dX2,dY2,dZ2,dX3,dY3,dZ3,dLambda,dRI,dVel;
double Rup,Rdown,FN,FM,refWeight;
double pdNCt[6];
Coords vdVel; /* Auxiliary velocity vector */
Coords eCRD_3DU2; /* Unit vector */
Coords eCRD_3DU3; /* Unit vector */
double dAbsorbed=0.0;
double pdGA[]={dGAy,dGAz};
double pdC[]={0.0,dDThetaHalf,dTlrnce,(double)abs_out+0.5};
coords_get(eCRD_3DU1,&dX1,&dY1,&dZ1); /* The components of the unit vector eCRD_3DU1 are constant */
pdTmp[10]=bb; /* The value of pdTmp[10] is constant throughout the run */
pdNCt[0]=x;
pdNCt[1]=y;
pdNCt[3]=vx;
pdNCt[4]=vy;
pdNCt[5]=vz;
if(lin<0.0) { /* Proximal polariser */
/*
First surface of the proximal polariser
=======================================
*/
switch(1) { /* A work-around enabling the propagation of neutrons which would otherwise be absorbed */
case 1:
pdNCt[2]=z-dZmax; /* Transformation accounting for the shift between the input and the own CSs */
pdC[0]=0.0;
DetermineIntersectionWithProximalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the first surface of the proximal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)ABSORB;
if(nType == 2)break;
/* A well-defined intersection with the first surface of the polariser has been identified; the question whether this intersection corresponds to an appropriate x value remains */
PROP_DT(dT);
/* Verify that the neutron hits the first surface of the polariser */
if(!IntersectsProximalPolariser(x,y,z-dZmax,h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
if(abs_out) { abs_reason[2]+=1; ABSORB; } else { break; }
} /* Endif(!IntersectsProximalPolariser(x,y,z-dZmax,h1,h2,dZmin,dZmax)) */
SCATTER; /* A well-defined intersection within the polariser bounds has been identified */
dSlp=-(bb*sin(dRt)+cos(dRt))/(bb*cos(dRt)-sin(dRt));
eCRD_3DU2=coords_set(0.0,dSlp/sqrt(1.0+dSlp*dSlp),1.0/sqrt(1.0+dSlp*dSlp));
eCRD_3DU3=coords_xp(eCRD_3DU2,eCRD_3DU1); /* Unit vector perpendicular to the surface of the proximal polariser (pointing towards the 'convex region') */
coords_get(eCRD_3DU2,&dX2,&dY2,&dZ2);
coords_get(eCRD_3DU3,&dX3,&dY3,&dZ3);
vdVel=coords_set(vx,vy,vz);
ObtainMagnitudeOfVector(vdVel,&dVel);
dTmp1=coords_sp(vdVel,eCRD_3DU1); /* Component of the velocity parallel to unit vector eCRD_3DU1 */
dTmp2=coords_sp(vdVel,eCRD_3DU2); /* Component of the velocity parallel to unit vector eCRD_3DU2 */
dTmp=dVel*dVel-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity */
/* Determine whether the neutron is reflected or refracted */
nRefracted=0;
isPolarising=0;
Rup=ObtainReflectivity(nRflctUp,dTmp3,R0,qc_Ni,alpha,m_u,W,m_residual,pTableUp);
Rdown=ObtainReflectivity(nRflctDn,dTmp3,R0,qc_Ni,alpha,m_d,W,m_residual,pTableDn);
if(Rup != Rdown) {
isPolarising = 1;
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
GetMonoPolRefProb(FN, FM, sx, &refWeight);
} else {
refWeight = Rup;
} /* Endif(Rup != Rdown) */
/* Verify that the weight refWeight is meaningful */
if(refWeight<0.0) { abs_reason[7] += 1; ABSORB; }
if(refWeight>1.0)refWeight=1.0;
if(rand01()>refWeight) { /* The neutron is refracted */
nRefracted=1;
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
} else { /* The neutron is reflected */
if (isPolarising) SetMonoPolRefOut(FN, FM, refWeight, &sy, &sx, &sz);
} /* Endif(rand01()>refWeight) */
/* Reflected neutrons are absorbed when abs_ref; the refraction is disabled when m_substrate == 0 or d_substrate == 0 */
if((!nRefracted) && abs_ref) { abs_reason[3] += 1; ABSORB; }
if(nRefracted && ((m_substrate<dTlrnce) || (d_substrate<dTlrnce)))break;
/* Proceed with the reflection or the refraction */
if(!nRefracted) { /* Reflection */
vx=dTmp1*dX1+dTmp2*dX2-dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2-dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp3*dZ3;
break;
} else { /* Refraction */
dLambda=2.0e10*PI*HBAR/MNEUTRON/dVel; /* [A] The de Broglie wavelength of the neutron */
dRI=(enable_ref) ? ObtainRefractiveIndex(dLambda,qc_substrate) : 1.0;
/* Obtain the velocity vector for the refracted neutron at the intersection */
dTmp=dRI*dRI*dVel*dVel-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
if(dTmp<0.0) { /* Total reflection */
if(abs_ref) { abs_reason[3] += 1; ABSORB; } else {
vx=dTmp1*dX1+dTmp2*dX2-dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2-dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp3*dZ3;
break;
} /* Endif(abs_ref) */
} /* Endif(dTmp<0.0) */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity of the refracted neutron */
/* The velocity vector for the refracted neutron at the intersection */
vx=dTmp1*dX1+dTmp2*dX2+dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2+dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp3*dZ3;
/*
Second surface of the proximal polariser
========================================
*/
pdNCt[0]=x;
pdNCt[1]=y;
pdNCt[2]=z-dZmax;
pdNCt[3]=vx;
pdNCt[4]=vy;
pdNCt[5]=vz;
pdC[0]=dOffsetY;
DetermineIntersectionWithProximalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the second surface of the proximal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)if(abs_single) { ABSORB; } else { break; }
if(nType == 2)break;
/* A well-defined intersection with the second surface of the polariser has been identified; the question whether this intersection corresponds to an appropriate x value remains */
PROP_DT(dT);
/* Verify that the neutron hits the second surface of the polariser */
if(!IntersectsProximalPolariser(x,y,z-dZmax,h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
if(abs_out) { abs_reason[6] += 1; ABSORB; } else { break; }
} /* Endif(!IntersectsProximalPolariser(x,y,z-dZmax,h1,h2,dZmin,dZmax)) */
SCATTER; /* A well-defined intersection within the polariser bounds has been identified */
dSlp=-(bb*sin(dRt)+cos(dRt))/(bb*cos(dRt)-sin(dRt));
eCRD_3DU2=coords_set(0.0,dSlp/sqrt(1.0+dSlp*dSlp),1.0/sqrt(1.0+dSlp*dSlp));
eCRD_3DU3=coords_xp(eCRD_3DU2,eCRD_3DU1); /* Unit vector perpendicular to the surface of the proximal polariser (pointing towards the 'convex region') */
coords_get(eCRD_3DU2,&dX2,&dY2,&dZ2);
coords_get(eCRD_3DU3,&dX3,&dY3,&dZ3);
vdVel=coords_set(vx,vy,vz);
ObtainMagnitudeOfVector(vdVel,&dVel);
dAbsorbed=1.0-exp(-dT*T_loss); // transmission losses due to absorption
dTmp1=coords_sp(vdVel,eCRD_3DU1); /* Component of the velocity parallel to unit vector eCRD_3DU1 */
dTmp2=coords_sp(vdVel,eCRD_3DU2); /* Component of the velocity parallel to unit vector eCRD_3DU2 */
dTmp=dVel*dVel-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity */
/* Special treatment in case that both surfaces are coated */
if(both_coated && (m_u>m_substrate) && (m_d>m_substrate)) {
isPolarising=0;
//Rup=ObtainReflectivity(nRflctUp,dTmp3,R0,qc_Ni,alpha,m_u-m_substrate,W,pTableUp);
//Rdown=ObtainReflectivity(nRflctDn,dTmp3,R0,qc_Ni,alpha,m_d-m_substrate,W,pTableDn);
if(Rup != Rdown) {
isPolarising = 1;
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
GetMonoPolRefProb(FN, FM, sx, &refWeight);
} else {
refWeight = Rup;
} /* Endif(Rup != Rdown) */
/* Verify that the weight refWeight is meaningful */
if(refWeight<0.0) { abs_reason[7] += 1; ABSORB; }
if(refWeight>1.0)refWeight=1.0;
if(rand01()>refWeight) { /* The neutron is refracted */
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
} else {
if(abs_ref) { abs_reason[7] += 1; ABSORB; } else {
/* The neutron is reflected, not yet implemented and should be done including multiple reflections */
vx=dTmp1*dX1+dTmp2*dX2-dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2-dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp3*dZ3;
break;
} /* Endif(abs_ref) */
} /* Endif(rand01()>refWeight) */
} /* Endif(both_coated) */
dLambda=2.0e10*PI*HBAR/MNEUTRON/dVel; /* [A] The de Broglie wavelength of the neutron */
dRI=(enable_ref) ? ObtainRefractiveIndex(dLambda,qc_substrate) : 1.0;
if(dRI == 0.0) { abs_reason[7] += 1; ABSORB; }
/* Obtain the velocity vector for the refracted neutron at the intersection */
dTmp=dVel*dVel/(dRI*dRI)-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity of the refracted neutron */
/* The velocity vector for the refracted neutron at the intersection */
vx=dTmp1*dX1+dTmp2*dX2+dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2+dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp3*dZ3;
} /* Endif(!nRefracted) */
} /* Endswitch(1) */
} else { /* Distal polariser */
/*
First surface of the distal polariser
=====================================
*/
switch(1) { /* A work-around enabling the propagation of neutrons which would otherwise be absorbed */
case 1:
pdNCt[2]=z+dZmin; /* Transformation accounting for the shift between the input and the own CSs */
pdC[0]=dOffsetY;
DetermineIntersectionWithDistalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the first surface of the distal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)ABSORB;
if(nType == 2)break;
/* A well-defined intersection with the first surface of the polariser has been identified; the question whether this intersection corresponds to an appropriate x value remains */
PROP_DT(dT);
/* Verify that the neutron hits the first surface of the polariser */
if(!IntersectsDistalPolariser(x,y,z+dZmin,h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
if(abs_out) { abs_reason[2]+=1; ABSORB; } else { break; }
} /* Endif(!IntersectsDistalPolariser(x,y,z+dZmin,h1,h2,dZmin,dZmax)) */
SCATTER; /* A well-defined intersection within the polariser bounds has been identified */
dSlp=(bb*sin(dRt)+cos(dRt))/(bb*cos(dRt)-sin(dRt));
eCRD_3DU2=coords_set(0.0,dSlp/sqrt(1.0+dSlp*dSlp),1.0/sqrt(1.0+dSlp*dSlp));
eCRD_3DU3=coords_xp(eCRD_3DU1,eCRD_3DU2); /* Unit vector perpendicular to the surface of the distal polariser (pointing towards the 'concave region') */
coords_get(eCRD_3DU2,&dX2,&dY2,&dZ2);
coords_get(eCRD_3DU3,&dX3,&dY3,&dZ3);
vdVel=coords_set(vx,vy,vz);
ObtainMagnitudeOfVector(vdVel,&dVel);
dTmp1=coords_sp(vdVel,eCRD_3DU1); /* Component of the velocity parallel to unit vector eCRD_3DU1 */
dTmp2=coords_sp(vdVel,eCRD_3DU2); /* Component of the velocity parallel to unit vector eCRD_3DU2 */
dTmp=dVel*dVel-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity */
/* Determine whether the neutron is reflected or refracted */
nRefracted=0;
isPolarising=0;
Rup=ObtainReflectivity(nRflctUp,dTmp3,R0,qc_Ni,alpha,m_u,W,m_residual,pTableUp);
Rdown=ObtainReflectivity(nRflctDn,dTmp3,R0,qc_Ni,alpha,m_d,W,m_residual,pTableDn);
if(Rup != Rdown) {
isPolarising = 1;
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
GetMonoPolRefProb(FN, FM, sx, &refWeight);
} else {
refWeight = Rup;
} /* Endif(Rup != Rdown) */
/* Verify that the weight refWeight is meaningful */
if(refWeight<0.0) { abs_reason[7] += 1; ABSORB; }
if(refWeight>1.0)refWeight=1.0;
if(rand01()>refWeight) { /* The neutron is refracted */
nRefracted=1;
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
} else { /* The neutron is reflected */
if (isPolarising) SetMonoPolRefOut(FN, FM, refWeight, &sy, &sx, &sz);
} /* Endif(rand01()>refWeight) */
/* Reflected neutrons are absorbed when abs_ref; the refraction is disabled when m_substrate == 0 or d_substrate == 0 */
if((!nRefracted) && abs_ref) { abs_reason[3] += 1; ABSORB; }
if(nRefracted && ((m_substrate<dTlrnce) || (d_substrate<dTlrnce)))break;
/* Proceed with the reflection or the refraction */
if(!nRefracted) { /* Reflection */
vx=dTmp1*dX1+dTmp2*dX2-dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2-dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp3*dZ3;
break;
} else { /* Refraction */
dLambda=2.0e10*PI*HBAR/MNEUTRON/dVel; /* [A] The de Broglie wavelength of the neutron */
dRI=(enable_ref) ? ObtainRefractiveIndex(dLambda,qc_substrate) : 1.0;
/* Obtain the velocity vector for the refracted neutron at the intersection */
dTmp=dRI*dRI*dVel*dVel-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
if(dTmp<0.0) { /* Total reflection */
if(abs_ref) { abs_reason[3] += 1; ABSORB; } else {
vx=dTmp1*dX1+dTmp2*dX2-dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2-dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp3*dZ3;
break;
} /* Endif(abs_ref) */
} /* Endif(dTmp<0.0) */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity of the refracted neutron */
/* The velocity vector for the refracted neutron at the intersection */
vx=dTmp1*dX1+dTmp2*dX2+dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2+dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp3*dZ3;
/*
Second surface of the distal polariser
======================================
*/
pdNCt[0]=x;
pdNCt[1]=y;
pdNCt[2]=z+dZmin;
pdNCt[3]=vx;
pdNCt[4]=vy;
pdNCt[5]=vz;
pdC[0]=0.0;
DetermineIntersectionWithDistalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the second surface of the distal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)if(abs_single) { ABSORB; } else { break; }
if(nType == 2)break;
/* A well-defined intersection with the second surface of the polariser has been identified; the question whether this intersection corresponds to an appropriate x value remains */
PROP_DT(dT);
/* Verify that the neutron hits the second surface of the polariser */
if(!IntersectsDistalPolariser(x,y,z+dZmin,h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
if(abs_out) { abs_reason[6] += 1; ABSORB; } else { break; }
} /* Endif(!IntersectsDistalPolariser(x,y,z+dZmin,h1,h2,dZmin,dZmax)) */
SCATTER; /* A well-defined intersection within the polariser bounds has been identified */
dSlp=(bb*sin(dRt)+cos(dRt))/(bb*cos(dRt)-sin(dRt));
eCRD_3DU2=coords_set(0.0,dSlp/sqrt(1.0+dSlp*dSlp),1.0/sqrt(1.0+dSlp*dSlp));
eCRD_3DU3=coords_xp(eCRD_3DU1,eCRD_3DU2); /* Unit vector perpendicular to the surface of the distal polariser (pointing towards the 'concave region') */
coords_get(eCRD_3DU2,&dX2,&dY2,&dZ2);
coords_get(eCRD_3DU3,&dX3,&dY3,&dZ3);
vdVel=coords_set(vx,vy,vz);
ObtainMagnitudeOfVector(vdVel,&dVel);
dAbsorbed=1.0-exp(-dT*T_loss); // transmission losses due to absorption
dTmp1=coords_sp(vdVel,eCRD_3DU1); /* Component of the velocity parallel to unit vector eCRD_3DU1 */
dTmp2=coords_sp(vdVel,eCRD_3DU2); /* Component of the velocity parallel to unit vector eCRD_3DU2 */
dTmp=dVel*dVel-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity; dTmp3 will be used later on */
/* Special treatment in case that both surfaces are coated */
if(both_coated && (m_u>m_substrate) && (m_d>m_substrate)) {
isPolarising=0;
//Rup=ObtainReflectivity(nRflctUp,dTmp3,R0,qc_Ni,alpha,m_u-m_substrate,W,pTableUp);
//Rdown=ObtainReflectivity(nRflctDn,dTmp3,R0,qc_Ni,alpha,m_d-m_substrate,W,pTableDn);
if(Rup != Rdown) {
isPolarising = 1;
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
GetMonoPolRefProb(FN, FM, sx, &refWeight);
} else {
refWeight = Rup;
} /* Endif(Rup != Rdown) */
/* Verify that the weight refWeight is meaningful */
if(refWeight<0.0) { abs_reason[7] += 1; ABSORB; }
if(refWeight>1.0)refWeight=1.0;
if(rand01()>refWeight) { /* The neutron is refracted */
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
} else {
if(abs_ref) { abs_reason[7] += 1; ABSORB; } else {
/* The neutron is reflected, not yet implemented and should be done including multiple reflections */
vx=dTmp1*dX1+dTmp2*dX2-dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2-dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp3*dZ3;
break;
} /* Endif(abs_ref) */
} /* Endif(rand01()>refWeight) */
} /* Endif(both_coated) */
dLambda=2.0e10*PI*HBAR/MNEUTRON/dVel; /* [A] The de Broglie wavelength of the neutron */
dRI=(enable_ref) ? ObtainRefractiveIndex(dLambda,qc_substrate) : 1.0;
if(dRI == 0.0) { abs_reason[7] += 1; ABSORB; }
/* Obtain the velocity vector for the refracted neutron at the intersection */
dTmp=dVel*dVel/(dRI*dRI)-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
dTmp3=sqrt(dTmp); /* Transverse component of the velocity of the refracted neutron */
/* The velocity vector for the refracted neutron at the intersection */
vx=dTmp1*dX1+dTmp2*dX2+dTmp3*dX3;
vy=dTmp1*dY1+dTmp2*dY2+dTmp3*dY3;
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp3*dZ3;
} /* Endif(!nRefracted) */
} /* Endswitch(1) */
} /* Endif(lin<0.0) */
/* Adjust the weight for absorption effects */
p*=1.-dAbsorbed;
/* Propagation of the neutron to the end of the component */
pdTmp[0]=z-(dZmax-dZmin);
pdTmp[1]=vz;
pdTmp[2]=0.5*dGAz;
EstimateExitTime(&dT);
if(dT <= 0.0) { abs_reason[7] += 1; ABSORB; } /* Estimation failed */
PROP_DT(dT);
%}
/**********************************************************************************************************************
De-allocation of memory */
FINALLY
%{
if(pTableDn)free(pTableDn);
if(pTableUp)free(pTableUp);
pTableDn=NULL;
pTableUp=NULL;
if(debug_print) {
fprintf(stdout,"abs_reason surface 1: %i, %i, %i, %i\n",abs_reason[0],abs_reason[1],abs_reason[2],abs_reason[3]);
fprintf(stdout,"abs_reason surface 2: %i, %i, %i, %i\n",abs_reason[4],abs_reason[5],abs_reason[6],abs_reason[7]);
} /* Endif(debug_print) */
%}
/**********************************************************************************************************************
Display of the polariser */
MCDISPLAY
%{
/* Re-initialisation of the parameters in case several instances are used */
dZmin=(lin<0.0) ? -(lin+length) : lin;
dZmax=(lin<0.0) ? -lin : lin+length;
ab=sqrt(dZmax*dZmin)/cos(dDThetaHalf); /* [m] The parameter a of the spiral curve */
bb=log(dZmax/dZmin)/delta_theta; /* [rad^-1] The parameter b of the spiral curve */
dOffsetY=(dZmin*sqrt(bb*bb+1.0)*(exp(bb*delta_theta)-1.0)*d_substrate/bb/(dZmax-dZmin));
int nNoPnts=101; /* Number of points used in the approximation of the (z,y) spiral */
int nCnt;
double dThetaPos=-0.5*delta_theta;
double xN[nNoPnts];
double yN[nNoPnts];
double zN[nNoPnts];
double dStp=delta_theta/(double)(nNoPnts-1);
for(nCnt=0;nCnt<nNoPnts;++nCnt) {
zN[nCnt]=-ab*exp(bb*dThetaPos)*cos(dThetaPos);
if(lin>0.0)zN[nCnt]=-zN[nCnt];
xN[nCnt]=0.5*((h2-h1)*fabs(zN[nCnt])+dZmax*h1-dZmin*h2)/(dZmax-dZmin);
if(lin>0.0) {
zN[nCnt] -= dZmin;
} else {
zN[nCnt] += dZmax;
}
yN[nCnt]=ab*exp(bb*dThetaPos)*sin(dThetaPos);
dThetaPos += dStp;
}
/* Display the polariser with three 'horizontal' curves and three vertical lines */
/* The three 'horizontal' curves */
for(nCnt=0;nCnt<nNoPnts-1;++nCnt)line(0.0,yN[nCnt],zN[nCnt],0.0,yN[nCnt+1],zN[nCnt+1]);
for(nCnt=0;nCnt<nNoPnts-1;++nCnt)line(xN[nCnt],yN[nCnt],zN[nCnt],xN[nCnt+1],yN[nCnt+1],zN[nCnt+1]);
for(nCnt=0;nCnt<nNoPnts-1;++nCnt)line(-xN[nCnt],yN[nCnt],zN[nCnt],-xN[nCnt+1],yN[nCnt+1],zN[nCnt+1]);
/* The three vertical curves */
line(xN[0],yN[0],zN[0],-xN[0],yN[0],zN[0]);
line(xN[(nNoPnts-1)/2],0.0,zN[(nNoPnts-1)/2],-xN[(nNoPnts-1)/2],0.0,zN[(nNoPnts-1)/2]);
line(xN[nNoPnts-1],yN[nNoPnts-1],zN[nNoPnts-1],-xN[nNoPnts-1],yN[nNoPnts-1],zN[nNoPnts-1]);
/* The three 'horizontal' curves */
for(nCnt=0;nCnt<nNoPnts-1;++nCnt)line(0.0,yN[nCnt]+dOffsetY,zN[nCnt],0.0,yN[nCnt+1]+dOffsetY,zN[nCnt+1]);
for(nCnt=0;nCnt<nNoPnts-1;++nCnt)line(xN[nCnt],yN[nCnt]+dOffsetY,zN[nCnt],xN[nCnt+1],yN[nCnt+1]+dOffsetY,zN[nCnt+1]);
for(nCnt=0;nCnt<nNoPnts-1;++nCnt)line(-xN[nCnt],yN[nCnt]+dOffsetY,zN[nCnt],-xN[nCnt+1],yN[nCnt+1]+dOffsetY,zN[nCnt+1]);
/* The three vertical curves */
line(xN[0],yN[0]+dOffsetY,zN[0],-xN[0],yN[0]+dOffsetY,zN[0]);
line(xN[(nNoPnts-1)/2],+dOffsetY,zN[(nNoPnts-1)/2],-xN[(nNoPnts-1)/2],+dOffsetY,zN[(nNoPnts-1)/2]);
line(xN[nNoPnts-1],yN[nNoPnts-1]+dOffsetY,zN[nNoPnts-1],-xN[nNoPnts-1],yN[nNoPnts-1]+dOffsetY,zN[nNoPnts-1]);
/* Start/End */
line(xN[0],yN[0],zN[0],xN[0],yN[0]+dOffsetY,zN[0]);
line(-xN[0],yN[0],zN[0],-xN[0],yN[0]+dOffsetY,zN[0]);
line(xN[nNoPnts-1],yN[nNoPnts-1],zN[nNoPnts-1],xN[nNoPnts-1],yN[nNoPnts-1]+dOffsetY,zN[nNoPnts-1]);
line(-xN[nNoPnts-1],yN[nNoPnts-1],zN[nNoPnts-1],-xN[nNoPnts-1],yN[nNoPnts-1]+dOffsetY,zN[nNoPnts-1]);
%}
END
Reference in New Issue View Git Blame Copy Permalink