1146 lines
52 KiB
Plaintext
1146 lines
52 KiB
Plaintext
/**********************************************************************************************************************
|
|
* McStas instrument definitions: http://www.mcstas.org
|
|
*
|
|
* Instrument: ESS Reflectometer Estia (ESS_Reflectometer_Estia)
|
|
* Component: Polariser
|
|
* Author: Evangelos Matsinos (evangelos.matsinos@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
|
|
* ----------------
|
|
* nIncRefr: [] 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)
|
|
* 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 nIncRefr = 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,
|
|
W = 0.0033333,
|
|
T_loss=0.125,
|
|
string reflect_d = "ReflectivitySpinDown.dat",
|
|
string reflect_u = "ReflectivitySpinUp.dat")
|
|
|
|
OUTPUT PARAMETERS (dDThetaHalf, dZmin, dZmax, ab, bb, dOffsetY, qc_substrate, nRflctUp, nRflctDn, localG, dGAx, dGAy, dGAz, dTmp)
|
|
|
|
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 nRflctDn,nRflctUp;
|
|
// double dTlrnce,dSpdLght,dMassNeut,dDThetaHalf,dhbarc,qc_Ni,qc_substrate,ab,bb,dOffsetY,dZmin,dZmax,dGAx,dGAy,dGAz;
|
|
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 abs_reason[]={0,0,0,0,0,0,0,0,0,0,0,0,0};
|
|
// Coords localG; /* Gravitational-acceleration vector */
|
|
// Coords eCRD_3DU1; /* Unit vector */
|
|
|
|
int debug_print=0;
|
|
|
|
/*
|
|
-----------------------------------------------------------------------------------------------------------------------
|
|
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);};
|
|
|
|
/*
|
|
-----------------------------------------------------------------------------------------------------------------------
|
|
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 positivity of the time interval between the initial state and the neutron incidence on the polariser.
|
|
Assertion of the appropriate reflection off the spiral surface of the proximal polariser.
|
|
*/
|
|
int IntersectsFirstPolariser(const double dT,
|
|
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(dT <= 0.0)return 0;
|
|
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 positivity of the time interval between the initial state and the neutron incidence on the polariser.
|
|
Assertion of the appropriate reflection off the spiral surface of the distal polariser.
|
|
*/
|
|
int IntersectsSecondPolariser(const double dT,
|
|
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(dT <= 0.0)return 0;
|
|
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 dMslope,
|
|
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)-(dMslope-1),-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;
|
|
|
|
}
|
|
|
|
%}
|
|
|
|
|
|
DECLARE
|
|
%{
|
|
|
|
int nRflctDn,nRflctUp;
|
|
double dTlrnce,dSpdLght,dMassNeut,dDThetaHalf,dhbarc,qc_Ni,qc_substrate,ab,bb,dOffsetY,dZmin,dZmax,dGAx,dGAy,dGAz;
|
|
Coords localG; /* Gravitational-acceleration vector */
|
|
Coords eCRD_3DU1; /* Unit vector */
|
|
|
|
int nSuccess,nRefracted,isPolarising;
|
|
double dRt,dT,dSlp,dWeight,dTmp,dTmp1,dTmp2,dTmpT,dX1,dY1,dZ1,dX2,dY2,dZ2,dX,dY,dZ,dLambda,dRI,dRnd,dVel;
|
|
double Rup,Rdown,FN,FM,refWeight;
|
|
Coords vdPos; /* Auxiliary position vector */
|
|
Coords vdVel; /* Auxiliary velocity vector */
|
|
Coords eCRD_3DU2; /* Unit vector */
|
|
Coords eCRD_3DU3; /* Unit vector */
|
|
|
|
%}
|
|
|
|
|
|
/**********************************************************************************************************************
|
|
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=(dZmin*sqrt(bb*bb+1.0)*(exp(bb*delta_theta)-1.0)*d_substrate/bb/(dZmax-dZmin)); /* [m] The offset of the inner
|
|
surface of the polariser in the y direction (with respect to the outer surface) */
|
|
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 along the x-axis */
|
|
|
|
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",nIncRefr);
|
|
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,"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) */
|
|
}
|
|
%}
|
|
|
|
|
|
/**********************************************************************************************************************
|
|
The main body of the method */
|
|
TRACE
|
|
%{
|
|
|
|
if(lin<0.0) { /* Proximal polariser */
|
|
|
|
// propagate to component entrance, will trace backward in case of two subsequent polarizers
|
|
//ALLOW_BACKPROP;
|
|
//PROP_Z0;
|
|
|
|
|
|
/*
|
|
Outer surface of the proximal polariser (towards the incoming beam)
|
|
===================================================================
|
|
*/
|
|
switch(1)
|
|
{
|
|
case 1:
|
|
|
|
if(dGAy == 0.0) { /* Case: no gravity along the y axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAz*ab*ab;
|
|
pdTmp[1]=ab*(vy*vz-dGAz*y);
|
|
pdTmp[2]=ab*vy*vy;
|
|
pdTmp[3]=vy*vy*(z-dZmax)-vy*vz*y+0.5*dGAz*y*y;
|
|
if(FTN1(-dDThetaHalf)*FTN1(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_out) { abs_reason[0]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN1,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[1]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(ab*exp(bb*dRt)*sin(dRt)-y)/vy; /* 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 along the z axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAy*ab*ab;
|
|
pdTmp[1]=ab*(dGAy*(z-dZmax)-vy*vz);
|
|
pdTmp[2]=-ab*vz*vz;
|
|
pdTmp[3]=vz*vz*y-vy*vz*(z-dZmax)+0.5*dGAy*(z-dZmax)*(z-dZmax);
|
|
if(FTN2(-dDThetaHalf)*FTN2(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_out) { abs_reason[0]+=1; ABSORB; } else { break;}
|
|
}
|
|
dRt=FindRoot(&FTN2,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[1]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=-(ab*exp(bb*dRt)*cos(dRt)+(z-dZmax))/vz; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
|
|
} else { /* Case: gravitational acceleration components along both the y and z axes (general case) */
|
|
pdTmp[0]=ab*dGAz;
|
|
pdTmp[1]=ab*dGAy;
|
|
pdTmp[2]=dGAy*(z-dZmax)-dGAz*y;
|
|
pdTmp[3]=dGAz*vy-dGAy*vz;
|
|
pdTmp[4]=y;
|
|
pdTmp[5]=vy;
|
|
pdTmp[6]=0.5*dGAy;
|
|
pdTmp[7]=(z-dZmax);
|
|
pdTmp[8]=vz;
|
|
pdTmp[9]=0.5*dGAz;
|
|
if(FTN3(-dDThetaHalf)*FTN3(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_out) { abs_reason[0]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN3,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[1]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(exp(bb*dRt)*(pdTmp[0]*sin(dRt)+pdTmp[1]*cos(dRt))+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) */
|
|
/* If the neutron continued propagating, an intersection with the mirror has been identified; the question whether this intersection corresponds to the surface of the polariser remains */
|
|
PROP_DT(dT);
|
|
vdVel=coords_set(vx,vy,vz); /* Velocity vector at the intersection */
|
|
/* Find whether the neutron hits the outer surface of the polariser */
|
|
if(!IntersectsFirstPolariser(dT,x,y,(z-dZmax),h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
|
|
if(abs_out) { abs_reason[2]+=1; ABSORB; } else { break; }
|
|
}
|
|
SCATTER;
|
|
|
|
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,&dX,&dY,&dZ);
|
|
if(!ObtainMagnitudeOfVector(vdVel,&dTmp)) { abs_reason[3]+=1; ABSORB; } /* Magnitude of the velocity */
|
|
dVel=dTmp; /* For use in case of refraction */
|
|
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=dTmp*dTmp-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
|
|
if(dTmp<0.0) { abs_reason[4]+=1; ABSORB; } /* Unphysical */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
dTmpT=dTmp; /* For the total reflection case */
|
|
|
|
/* Determine whether the neutron is reflected or refracted */
|
|
nRefracted=0;
|
|
isPolarising=0;
|
|
|
|
Rup=ObtainReflectivity(nRflctUp,dTmp,R0,qc_Ni,alpha,m_u,W,m_substrate,pTableUp);
|
|
Rdown=ObtainReflectivity(nRflctDn,dTmp,R0,qc_Ni,alpha,m_d,W,m_substrate,pTableDn);
|
|
if(Rup != Rdown) {
|
|
isPolarising = 1;
|
|
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
|
|
GetMonoPolRefProb(FN, FM, sx, &refWeight);
|
|
} else
|
|
refWeight = Rup;
|
|
|
|
// check that refWeight is meaningfull
|
|
if (refWeight < 0) ABSORB;
|
|
if (refWeight > 1) refWeight =1 ;
|
|
|
|
if(rand01()>refWeight) {
|
|
// neutron is transmitted
|
|
nRefracted=1;
|
|
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
} else {
|
|
// neutron is reflected
|
|
if (isPolarising) SetMonoPolRefOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
}
|
|
|
|
|
|
/* Reflected neutrons are absorbed when abs_ref; the refraction is disabled when m_substrate == 0 */
|
|
if((!nRefracted) && abs_ref) { abs_reason[5]+=1; ABSORB; }
|
|
if(nRefracted && (m_substrate<dTlrnce)) { break; }
|
|
/* Proceed with the reflection or the refraction */
|
|
if(!nRefracted) { /* Reflection */
|
|
vx=dTmp1*dX1+dTmp2*dX2-dTmpT*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2-dTmpT*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2-dTmpT*dZ;
|
|
break;
|
|
} else { /* Refraction */
|
|
dLambda=2.0e-05*PI*dhbarc*dSpdLght/dMassNeut/dVel; /* [A] The de Broglie wavelength of the neutron */
|
|
dRI=(nIncRefr) ? 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[6]+=1; ABSORB;
|
|
} else {
|
|
vx=dTmp1*dX1+dTmp2*dX2-dTmpT*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2-dTmpT*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2-dTmpT*dZ;
|
|
break;
|
|
} /* Endif(abs_ref) */
|
|
} /* Endif(dTmp<0.0) */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
/* The velocity vector for the refracted neutron at the intersection */
|
|
vx=dTmp1*dX1+dTmp2*dX2+dTmp*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2+dTmp*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp*dZ;
|
|
} /* Endif(!nRefracted) */
|
|
}
|
|
/*
|
|
Inner surface of the proximal polariser (away from the incoming beam)
|
|
=====================================================================
|
|
*/
|
|
dT=(dZmax-z-dZmin)/vz; // time to propagate to end of component (if second surface is not reached
|
|
switch(nRefracted)
|
|
{
|
|
case 1:
|
|
p*=1.0-T_loss; // absorb some neutrons in the substrate
|
|
if(dGAy == 0.0) { /* Case: no gravity along the y axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAz*ab*ab;
|
|
pdTmp[1]=ab*(vy*vz-dGAz*(y-dOffsetY));
|
|
pdTmp[2]=ab*vy*vy;
|
|
pdTmp[3]=vy*vy*(z-dZmax)-vy*vz*(y-dOffsetY)+0.5*dGAz*(y-dOffsetY)*(y-dOffsetY);
|
|
if(FTN1(-dDThetaHalf)*FTN1(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
PROP_DT(dT);
|
|
if(abs_single) { abs_reason[7]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN1,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
PROP_DT(dT);
|
|
if(abs_out) { abs_reason[8]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(ab*exp(bb*dRt)*sin(dRt)-(y-dOffsetY))/vy; /* 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 along the z axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAy*ab*ab;
|
|
pdTmp[1]=ab*(dGAy*(z-dZmax)-vy*vz);
|
|
pdTmp[2]=-ab*vz*vz;
|
|
pdTmp[3]=vz*vz*(y-dOffsetY)-vy*vz*(z-dZmax)+0.5*dGAy*(z-dZmax)*(z-dZmax);
|
|
if(FTN2(-dDThetaHalf)*FTN2(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
PROP_DT(dT);
|
|
if(abs_single) { abs_reason[7]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN2,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
PROP_DT(dT);
|
|
if(abs_out) { abs_reason[8]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=-(ab*exp(bb*dRt)*cos(dRt)+(z-dZmax))/vz; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
|
|
} else { /* Case: gravitational acceleration components along both the y and z axes (general case) */
|
|
pdTmp[0]=ab*dGAz;
|
|
pdTmp[1]=ab*dGAy;
|
|
pdTmp[2]=dGAy*(z-dZmax)-dGAz*(y-dOffsetY);
|
|
pdTmp[3]=dGAz*vy-dGAy*vz;
|
|
pdTmp[4]=y-dOffsetY;
|
|
pdTmp[5]=vy;
|
|
pdTmp[6]=0.5*dGAy;
|
|
pdTmp[7]=(z-dZmax);
|
|
pdTmp[8]=vz;
|
|
pdTmp[9]=0.5*dGAz;
|
|
if(FTN3(-dDThetaHalf)*FTN3(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
PROP_DT(dT);
|
|
if(abs_single) { abs_reason[7]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN3,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
PROP_DT(dT);
|
|
if(abs_out) { abs_reason[8]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(exp(bb*dRt)*(pdTmp[0]*sin(dRt)+pdTmp[1]*cos(dRt))+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) */
|
|
/* If the neutron continued propagating, an intersection with the mirror has been identified; the question whether this intersection corresponds to the surface of the real polariser remains */
|
|
PROP_DT(dT);
|
|
vdVel=coords_set(vx,vy,vz); /* Velocity vector at the intersection */
|
|
/* Find whether the neutron hits the inner surface of the polariser */
|
|
if(!IntersectsFirstPolariser(dT,x,(y-dOffsetY),(z-dZmax),h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
|
|
if(abs_out) { abs_reason[9]+=1; ABSORB; } else { break; }
|
|
}
|
|
SCATTER;
|
|
|
|
if (both_coated) {
|
|
nRefracted=0;
|
|
isPolarising=0;
|
|
|
|
Rup=ObtainReflectivity(nRflctUp,dTmpT,R0,qc_Ni,alpha,m_u,W,0.0,pTableUp);
|
|
Rdown=ObtainReflectivity(nRflctDn,dTmpT,R0,qc_Ni,alpha,m_d,W,0.0,pTableDn);
|
|
if(Rup != Rdown) {
|
|
isPolarising = 1;
|
|
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
|
|
GetMonoPolRefProb(FN, FM, sx, &refWeight);
|
|
} else
|
|
refWeight = Rup;
|
|
|
|
// check that refWeight is meaningfull
|
|
if (refWeight < 0) ABSORB;
|
|
if (refWeight > 1) refWeight =1 ;
|
|
|
|
if(rand01()>refWeight) {
|
|
// neutron is transmitted
|
|
nRefracted=1;
|
|
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
} else {
|
|
ABSORB;
|
|
// neutron is reflected, not yet implemented and should be done including multiple reflections
|
|
if (isPolarising) SetMonoPolRefOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
if(!ObtainMagnitudeOfVector(vdVel,&dTmp)){ abs_reason[10]+=1; ABSORB; } /* Magnitude of the velocity */
|
|
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=dTmp*dTmp-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
|
|
if(dTmp<0.0) { abs_reason[12]+=1; ABSORB; } /* Unphysical */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
vx=dTmp1*dX1+dTmp2*dX2-dTmp*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2-dTmp*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp*dZ;
|
|
break;
|
|
}
|
|
}
|
|
|
|
/*
|
|
Propagation of the neutron to the air
|
|
=====================================
|
|
*/
|
|
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') */
|
|
if(!ObtainMagnitudeOfVector(vdVel,&dTmp)){ abs_reason[10]+=1; ABSORB; } /* Magnitude of the velocity */
|
|
dLambda=2.0e-05*PI*dhbarc*dSpdLght/dMassNeut/dTmp; /* [A] The de Broglie wavelength of the neutron */
|
|
dRI=(nIncRefr) ? ObtainRefractiveIndex(dLambda,qc_substrate) : 1.0;
|
|
if(dRI<dTlrnce) { abs_reason[11]+=1; ABSORB; } /* Unphysical */
|
|
/* Obtain the velocity vector for the refracted neutron at the intersection */
|
|
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=dTmp*dTmp/(dRI*dRI)-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
|
|
if(dTmp<0.0) { abs_reason[12]+=1; ABSORB; } /* Unphysical */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
/* The velocity vector for the refracted neutron at the intersection */
|
|
coords_get(eCRD_3DU2,&dX2,&dY2,&dZ2);
|
|
coords_get(eCRD_3DU3,&dX,&dY,&dZ);
|
|
vx=dTmp1*dX1+dTmp2*dX2+dTmp*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2+dTmp*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp*dZ;
|
|
}
|
|
|
|
} else { /* Distal polariser */
|
|
|
|
/*
|
|
Inner surface of the distal polariser (towards the incoming beam)
|
|
=================================================================
|
|
*/
|
|
switch(1)
|
|
{
|
|
case 1:
|
|
if(dGAy == 0.0) { /* Case: no gravity along the y axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAz*ab*ab;
|
|
pdTmp[1]=ab*(vy*vz-dGAz*(y-dOffsetY));
|
|
pdTmp[2]=-ab*vy*vy;
|
|
pdTmp[3]=vy*vy*(z+dZmin)-vy*vz*(y-dOffsetY)+0.5*dGAz*(y-dOffsetY)*(y-dOffsetY);
|
|
if(FTN1(-dDThetaHalf)*FTN1(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_single) { abs_reason[0]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN1,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) {abs_reason[1]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(ab*exp(bb*dRt)*sin(dRt)-(y-dOffsetY))/vy; /* 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 along the z axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAy*ab*ab;
|
|
pdTmp[1]=ab*(vy*vz-dGAy*(z+dZmin));
|
|
pdTmp[2]=-ab*vz*vz;
|
|
pdTmp[3]=vz*vz*(y-dOffsetY)-vy*vz*(z+dZmin)+0.5*dGAy*(z+dZmin)*(z+dZmin);
|
|
if(FTN2(-dDThetaHalf)*FTN2(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_single) { abs_reason[0]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN2,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[1]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(ab*exp(bb*dRt)*cos(dRt)-z)/vz; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
|
|
} else { /* Case: gravitational acceleration components along both the y and z axes (general case) */
|
|
pdTmp[0]=ab*dGAz;
|
|
pdTmp[1]=-ab*dGAy;
|
|
pdTmp[2]=dGAy*(z+dZmin)-dGAz*(y-dOffsetY);
|
|
pdTmp[3]=dGAz*vy-dGAy*vz;
|
|
pdTmp[4]=y-dOffsetY;
|
|
pdTmp[5]=vy;
|
|
pdTmp[6]=0.5*dGAy;
|
|
pdTmp[7]=-(z+dZmin);
|
|
pdTmp[8]=-vz;
|
|
pdTmp[9]=-0.5*dGAz;
|
|
if(FTN3(-dDThetaHalf)*FTN3(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_single) { abs_reason[0]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN3,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[1]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(exp(bb*dRt)*(pdTmp[0]*sin(dRt)+pdTmp[1]*cos(dRt))+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) */
|
|
/* If the neutron continued propagating, an intersection with the mirror has been identified; the question whether this intersection corresponds to the surface of the real polariser remains */
|
|
PROP_DT(dT);
|
|
vdVel=coords_set(vx,vy,vz); /* Velocity vector at the intersection */
|
|
/* Find whether the neutron hits the inner surface of the polariser */
|
|
if(!IntersectsSecondPolariser(dT,x,(y-dOffsetY),(z+dZmin),h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
|
|
if(abs_out) { abs_reason[2]+=1; ABSORB; } else { break; }
|
|
}
|
|
SCATTER;
|
|
|
|
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,&dX,&dY,&dZ);
|
|
if(!ObtainMagnitudeOfVector(vdVel,&dTmp)) { abs_reason[3]+=1; ABSORB; } /* Magnitude of the velocity */
|
|
dVel=dTmp; /* For use in case of refraction */
|
|
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=dTmp*dTmp-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
|
|
if(dTmp<0.0) { abs_reason[4]+=1; ABSORB; } /* Unphysical */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
dTmpT=dTmp; /* For the total reflection case */
|
|
/* Determine whether the neutron is reflected or refracted */
|
|
nRefracted=0;
|
|
isPolarising=0;
|
|
|
|
Rup=ObtainReflectivity(nRflctUp,dTmp,R0,qc_Ni,alpha,m_u,W,m_substrate,pTableUp);
|
|
Rdown=ObtainReflectivity(nRflctDn,dTmp,R0,qc_Ni,alpha,m_d,W,m_substrate,pTableDn);
|
|
if(Rup != Rdown) {
|
|
isPolarising = 1;
|
|
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
|
|
GetMonoPolRefProb(FN, FM, sx, &refWeight);
|
|
} else
|
|
refWeight = Rup;
|
|
|
|
// check that refWeight is meaningfull
|
|
if (refWeight < 0) ABSORB;
|
|
if (refWeight > 1) refWeight =1 ;
|
|
|
|
if(rand01()>refWeight) {
|
|
// neutron is transmitted
|
|
nRefracted=1;
|
|
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
} else {
|
|
// neutron is reflected
|
|
if (isPolarising) SetMonoPolRefOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
}
|
|
|
|
/* Reflected neutrons are absorbed when abs_ref; the refraction is disabled when m_substrate == 0 */
|
|
if((!nRefracted) && abs_ref) { abs_reason[5]+=1; ABSORB; }
|
|
if(nRefracted && (m_substrate<dTlrnce)) { SCATTER; }
|
|
/* Proceed with the reflection or the refraction */
|
|
if(!nRefracted) { /* Reflection */
|
|
vx=dTmp1*dX1+dTmp2*dX2-dTmp*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2-dTmp*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2-dTmp*dZ;
|
|
SCATTER;
|
|
} else { /* Refraction */
|
|
dLambda=2.0e-05*PI*dhbarc*dSpdLght/dMassNeut/dVel; /* [A] The de Broglie wavelength of the neutron */
|
|
dRI=(nIncRefr) ? 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[6]+=1; ABSORB;
|
|
} else {
|
|
vx=dTmp1*dX1+dTmp2*dX2-dTmpT*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2-dTmpT*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2-dTmpT*dZ;
|
|
SCATTER;
|
|
} /* Endif(abs_ref) */
|
|
} /* Endif(dTmp<0.0) */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
/* The velocity vector for the refracted neutron at the intersection */
|
|
vx=dTmp1*dX1+dTmp2*dX2+dTmp*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2+dTmp*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp*dZ;
|
|
} /* Endif(!nRefracted) */
|
|
}
|
|
/*
|
|
Outer surface of the distal polariser (away from the incoming beam)
|
|
===================================================================
|
|
*/
|
|
switch(nRefracted)
|
|
{
|
|
case 1:
|
|
p*=1.0-T_loss; // absorb some neutrons in the substrate
|
|
if(dGAy == 0.0) { /* Case: no gravity along the y axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAz*ab*ab;
|
|
pdTmp[1]=ab*(vy*vz-dGAz*y);
|
|
pdTmp[2]=-ab*vy*vy;
|
|
pdTmp[3]=vy*vy*(z+dZmin)-vy*vz*y+0.5*dGAz*y*y;
|
|
if(FTN1(-dDThetaHalf)*FTN1(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_out) { abs_reason[7]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN1,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[8]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(ab*exp(bb*dRt)*sin(dRt)-y)/vy; /* 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 along the z axis (simplified case) */
|
|
pdTmp[0]=0.5*dGAy*ab*ab;
|
|
pdTmp[1]=ab*(vy*vz-dGAy*(z+dZmin));
|
|
pdTmp[2]=-ab*vz*vz;
|
|
pdTmp[3]=vz*vz*y-vy*vz*(z+dZmin)+0.5*dGAy*(z+dZmin)*(z+dZmin);
|
|
if(FTN2(-dDThetaHalf)*FTN2(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_out) { abs_reason[7]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN2,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[8]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(ab*exp(bb*dRt)*cos(dRt)-(z+dZmin))/vz; /* The time needed for the neutron to travel from the original point (x,y,z) to the intersection */
|
|
} else { /* Case: gravitational acceleration components along both the y and z axes (general case) */
|
|
pdTmp[0]=ab*dGAz;
|
|
pdTmp[1]=-ab*dGAy;
|
|
pdTmp[2]=dGAy*(z+dZmin)-dGAz*y;
|
|
pdTmp[3]=dGAz*vy-dGAy*vz;
|
|
pdTmp[4]=y;
|
|
pdTmp[5]=vy;
|
|
pdTmp[6]=0.5*dGAy;
|
|
pdTmp[7]=-(z+dZmin);
|
|
pdTmp[8]=-vz;
|
|
pdTmp[9]=-0.5*dGAz;
|
|
if(FTN3(-dDThetaHalf)*FTN3(dDThetaHalf)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
|
|
if(abs_out) { abs_reason[7]+=1; ABSORB; } else { break; }
|
|
}
|
|
dRt=FindRoot(&FTN3,-dDThetaHalf,dDThetaHalf,dTlrnce,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
|
|
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on the abs_out variable */
|
|
if(abs_out) { abs_reason[8]+=1; ABSORB; } else { break; }
|
|
}
|
|
dT=(exp(bb*dRt)*(pdTmp[0]*sin(dRt)+pdTmp[1]*cos(dRt))+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) */
|
|
/* If the neutron continued propagating, an intersection with the mirror has been identified; the question whether this intersection corresponds to the surface of the real polariser remains */
|
|
PROP_DT(dT);
|
|
vdVel=coords_set(vx,vy,vz); /* Velocity vector at the intersection */
|
|
/* Find whether the neutron hits the outer surface of the polariser */
|
|
if(!IntersectsSecondPolariser(dT,x,y,(z+dZmin),h1,h2,dZmin,dZmax)) { /* Evidently, x out of bounds */
|
|
if(abs_out) { abs_reason[9]+=1; ABSORB; } else { break; }
|
|
}
|
|
SCATTER;
|
|
|
|
if (both_coated) {
|
|
nRefracted=0;
|
|
isPolarising=0;
|
|
|
|
Rup=ObtainReflectivity(nRflctUp,dTmpT,R0,qc_Ni,alpha,m_u,W,0.0,pTableUp);
|
|
Rdown=ObtainReflectivity(nRflctDn,dTmpT,R0,qc_Ni,alpha,m_d,W,0.0,pTableDn);
|
|
if(Rup != Rdown) {
|
|
isPolarising = 1;
|
|
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
|
|
GetMonoPolRefProb(FN, FM, sx, &refWeight);
|
|
} else
|
|
refWeight = Rup;
|
|
|
|
// check that refWeight is meaningfull
|
|
if (refWeight < 0) ABSORB;
|
|
if (refWeight > 1) refWeight =1 ;
|
|
|
|
if(rand01()>refWeight) {
|
|
// neutron is transmitted
|
|
nRefracted=1;
|
|
if (isPolarising) SetMonoPolTransOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
} else {
|
|
// neutron is reflected, not yet implemented
|
|
ABSORB;
|
|
if (isPolarising) SetMonoPolRefOut(FN, FM, refWeight, &sy, &sx, &sz);
|
|
}
|
|
}
|
|
|
|
/*
|
|
Propagation of the neutron to the air
|
|
=====================================
|
|
*/
|
|
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') */
|
|
if(!ObtainMagnitudeOfVector(vdVel,&dTmp)) { abs_reason[10]+=1; ABSORB; } /* Magnitude of the velocity */
|
|
dLambda=2.0e-05*PI*dhbarc*dSpdLght/dMassNeut/dTmp; /* [A] The de Broglie wavelength of the neutron */
|
|
dRI=(nIncRefr) ? ObtainRefractiveIndex(dLambda,qc_substrate) : 1.0;
|
|
if(dRI<dTlrnce) { abs_reason[11]+=1; ABSORB; } /* Unphysical */
|
|
/* Obtain the velocity vector for the refracted neutron at the intersection */
|
|
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=dTmp*dTmp/(dRI*dRI)-dTmp1*dTmp1-dTmp2*dTmp2; /* Square of the transverse component of the velocity */
|
|
if(dTmp<0.0) { abs_reason[12]+=1; ABSORB; } /* Unphysical */
|
|
dTmp=sqrt(dTmp); /* Transverse component of the velocity */
|
|
/* The velocity vector for the refracted neutron at the intersection */
|
|
coords_get(eCRD_3DU2,&dX2,&dY2,&dZ2);
|
|
coords_get(eCRD_3DU3,&dX,&dY,&dZ);
|
|
vx=dTmp1*dX1+dTmp2*dX2+dTmp*dX;
|
|
vy=dTmp1*dY1+dTmp2*dY2+dTmp*dY;
|
|
vz=dTmp1*dZ1+dTmp2*dZ2+dTmp*dZ;
|
|
|
|
}
|
|
} /* Endif(lin<0.0) */
|
|
|
|
SCATTER;
|
|
|
|
%}
|
|
|
|
|
|
/**********************************************************************************************************************
|
|
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, %i, %i, %i\n",
|
|
abs_reason[0],abs_reason[1],abs_reason[2],abs_reason[3],abs_reason[4],abs_reason[5],abs_reason[6]);
|
|
fprintf(stdout,"abs_reason surface 2: %i,%i,%i,%i,%i,%i\n",
|
|
abs_reason[7],abs_reason[8],abs_reason[9],abs_reason[10],abs_reason[11],abs_reason[12]);
|
|
}
|
|
|
|
%}
|
|
|
|
|
|
/**********************************************************************************************************************
|
|
Display of the polariser */
|
|
MCDISPLAY
|
|
%{
|
|
|
|
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
|