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

211 lines
8.1 KiB
Plaintext
Raw Permalink Blame History

/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Scatter_log_iterator.comp
*
* %I
*
* Written by: Erik B Knudsen
* Updated by: Rodon Kolevatov
* Date: November 2012
* Version: $Revision: 1.21 $
* Release: McStas 2.1
* Origin: DTU Physics
*
* Iteration element for a Scatter_log
*
* %D
*
* This component marks the beginning of the region in trace in which pseudo-neutrons
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
* a set of Scatter_logger/Scatter_logger_stop components.
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
* as if they were regular neutrons in the classical McStas-manner.
*
* %P
* Input parameters:
*
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
*
* %E
*******************************************************************************/
DEFINE COMPONENT Shielding_log_iterator_Ni_new
DEFINITION PARAMETERS (compute_func=NULL)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (nstate_initial,s0,s1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
#ifndef M1THICKNESS
#define M1THICKNESS 1500.0 //thickness of the m=1 coating
#endif
/*This is the specialized pseudo-neutron function that computes
an escaping neutron from logged before and after SCATTER neutron states
with weight corresponding to capture in the Nickel layers*/
int exit_neutron_Ni(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
/*!!Note that the transformation into global coordinate system must be done while logging
as we do not have access to neither the component name nor can get the component rotation by index.*/
Coords c1,c2;
Rotation R1,R2;
/*so now compute the pseudo neutron state and possibly user variables*/
/*momentum transfer at reflection. ASSUMES NO GRAVITY???*/
double velocity=sqrt((S1->_vx)*(S1->_vx)+(S1->_vy)*(S1->_vy)+(S1->_vz)*(S1->_vz));
double qtransf = V2Q*sqrt((S1->_vx-S0->_vx)*(S1->_vx-S0->_vx)+(S1->_vy-S0->_vy)*(S1->_vy-S0->_vy)+(S1->_vz-S0->_vz)*(S1->_vz-S0->_vz));
double qinm = qtransf/0.0218;
/*position comes from "new" state*/
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
/*velocity is the "old" state*/
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
/*time from new*/
ns_tilde[6]=S1->_t;
/*spin comes from "new" state*/
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
/*weight of capture in Ni is determined analytically*/
double captureprob;
if ((S1->_mvalue<0)||(qinm<0.0001)) {ns_tilde[10]=0.0; captureprob=0.0;}
else if (qinm<=1.01){
/*q<qc(Ni), loss due to diffusion beyond the coating cutoff*/
/* captureprob = 4.0582E-09*2200.0/velocity/((0.085787L*18.5 + 0.008321L*5.7)*1.0E-08+4.0582E-09*2200.0/velocity)*(S0->_p-S1->_p)/S0->_p; */
double sigmadif ;//diffuse scattering cross section
if (velocity > 1950.0){ sigmadif = 18.5; //lambda < 2A, totally incoherent scattering
} else if (velocity < 1300.0){sigmadif = 5.2; // lambda >3A, totally coherent regime, only incoherent part contributes
} else sigmadif=5.2+13.3*(1950-velocity)/(1950.0-1300.0);
/*Now checking if the coating is m=1 or has higher m*/
if(S1->_mvalue<1.05) /*m=1 coating=> "triple thickness approximation" with Im k determined for reflectivity dip*/
{
double imk; /*imaginary part of momentum in the layer*/
imk=1.e-8*(sigmadif*0.09121+0.41*2200.0/velocity)*velocity/3956.0*M1THICKNESS; /*v(m/s)x\lambda(AA)=3956*/
captureprob = 4.49*2200.0/velocity/(sigmadif+4.49*2200.0/velocity)*(S0->_p-S1->_p)*(1-exp(-2.0*imk*M1THICKNESS*3.0))/S0->_p;
} else { // reflection loss below qcNi in case of multilayer.
/*Taking the minimum of what is for m=1 coating and what is when assume that neutron is physically lost beyond the cutoff.*/
double imk; /*imaginary part of momentum in the layer*/
imk=1.e-8*(sigmadif*0.09121+0.41*2200.0/velocity)*velocity/3956.0*M1THICKNESS; /*v(m/s)x\lambda(AA)=3956*/
double c1, c2;
c1 = 4.49*2200.0/velocity/(sigmadif+4.49*2200.0/velocity)*(S0->_p-S1->_p)*(1-exp(-2.0*imk*M1THICKNESS*3.0))/S0->_p;
c2 = 0.005*(S1->_mvalue+0.1)*(S0->_p-S1->_p)/S0->_p; // for the matters of conservative estimate might be divided by 1-R at mvalue+0.1.
captureprob = (c1>c2) ? c1 : c2 ;
}
} //end of "if (qinm<=1.02)"
else if (qinm<=S1->_mvalue+0.1) { /*q>q_c(Ni) and reflection, absorption per hit */
captureprob=0.005*qinm;
} else { /* transmission beyond smirrorcutoff, some reserve for m=1 coating*/
if(S1->_mvalue<1.05) captureprob=0.0025*1.3*1.3/qinm;
else captureprob=0.0025*(S1->_mvalue+0.1)*(S1->_mvalue+0.1)/qinm;}
/*check change in weight*/
if((S0->_p-S1->_p)>1.e-5*S0->_p) ns_tilde[10]=S0->_p*captureprob; else ns_tilde[10]=0.0;
/* if "mvalue" is -1, then absorption happened not on the coating but somewhere else. The weight of capture in nickel should be set to zero.*/
if(ns_tilde[10]>(S0->_p-S1->_p)) { printf("%f\t%f\t%f\t%E\t%E\t%f\t%E\n",velocity, qinm, S1->_mvalue, S0->_p, S1->_p, captureprob, ns_tilde[10]); exit(0);}
return 0;
}
#define NOABS \
do {/* Nothing*/} while(0)
%}
DECLARE
%{
int (*pseudo_neutron_state_function_Ni) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
struct Generalized_State_t *s1,*s0;
double *nstate_initial;
int optics_not_hit;
/*need a pointer to the structure set up by the logger*/
%}
INITIALIZE
%{
if (compute_func) {
pseudo_neutron_state_function_Ni=compute_func;
}else{
pseudo_neutron_state_function_Ni=exit_neutron_Ni;
}
nstate_initial=NULL;
optics_not_hit=0;
%}
TRACE
%{
//printf("Entering Iterator_Ni\n");
/*I am the start of the pseudo neutron iterator*/
if (nstate_initial==NULL){
optics_not_hit=0; /* Fresh start, resetting variable */
double *ns=nstate_initial=calloc(11,sizeof(double));
ns[0]=x;ns[1]=y; ns[2]=z;
ns[3]=vx;ns[4]=vy;ns[5]=vz;
ns[6]=t;
ns[7]=sx;ns[8]=sy;ns[9]=sz;
ns[10]=p;
s0=Bounce_store;
s1=Bounce_store+1;
/* Remove std. ABSORB to avoid breaking analysis loop */
#undef mcabsorb
#define mcabsorb scatter_iterator_stop
}
//if (s1->_p!=-1){
if (s1->_p!=-1 && s1->_p!=-2){
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
double nstate[11];
if ( pseudo_neutron_state_function_Ni(nstate,s0,s1) ){
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
}
/*set neutron state for subsequent components*/
x=nstate[0];y=nstate[1];z=nstate[2];
vx=nstate[3];vy=nstate[4];vz=nstate[5];
t=nstate[6];
sx=nstate[7];sy=nstate[8];sz=nstate[9];
p=nstate[10];
s0++;
s1++;
// printf("Ni: z=%g\tp=%g\n",z,p);
}else if (Bounce_store[1]._p==-1){
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
x=s0->_x;y=s0->_y;z=s0->_z;
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
t=s0->_t;
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
p=s0->_p;
optics_not_hit = 1;
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
} else {
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
fprintf(stderr,"This should not happen. Period.\n");
exit(1);
}
%}
MCDISPLAY
%{
/* A bit ugly; hard-coded dimensions. */
magnify("");
line(0,0,0,0.2,0,0);
line(0,0,0,0,0.2,0);
line(0,0,0,0,0,0.2);
%}
END
Reference in New Issue View Git Blame Copy Permalink