Estia-McStas/simulation/Shielding_log_iterator_total.comp

/*******************************************************************************
*
* 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 
* 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_total
  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 
%{
  /*This is the specialized pseudo-neutron function that computes
    an escaping neutron from logged before and after SCATTER neutron states*/
  int exit_neutron(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*/
    /*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 is difference old-new to mean the neutrons "deposited" in the guide wall*/
    ns_tilde[10]=S0->_p-S1->_p;
    return 0;
  }

#define NOABS						     \
  do {/* Nothing*/} while(0)
  
%}

DECLARE
%{
  int (*pseudo_neutron_state_function) (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=compute_func;
  }else{
    pseudo_neutron_state_function=exit_neutron;
  }
  nstate_initial=NULL;
    
  optics_not_hit=0;


%}


TRACE
%{
//printf("Entering Iterator_total\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(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("Tot: 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