Comp updates etc from ISIS

mcstas-comps/contrib/Commodus_I3.comp

@@ -0,0 +1,115 @@
+/*******************************************************************************
+* 
+*  McStas, neutron ray-tracing package
+*          Copyright (C) 1997-2008, All rights reserved
+*          Risoe National Laboratory, Roskilde, Denmark
+*          Institut Laue Langevin, Grenoble, France
+* 
+*  Component: Commodus_I3
+* 
+* 
+*  %I
+* Written by: G. Skoro, based on ViewModerator4 from S. Ansell
+*  Date: July 2022
+*  Origin: ISIS
+*
+* ISIS Moderators (Tested with McStas 3.1 (Windows))
+* 
+*  %D
+*  Produces a neutron distribution at the ISIS TS1 or TS2 corresponding moderator front face position.  
+*  The Face argument determines which TS1 or TS2 beamline is to be sampled by using corresponding file.
+*  Neutrons are created having a range of energies determined by the E0 and E1 arguments.
+*  Trajectories are produced such that they pass through the moderator face (defined by 
+*  modXsize and modZsize) and a focusing rectangle (defined by xw, yh and dist). 
+*  		--- HOW TO USE ---
+* 
+*       Example:   Commodus_I3(Face="TS1verBase2016_LH8020_newVM-var_South04_Merlin.mcstas", E0 = E_min, E1 = E_max,
+*       modXsize = 0.12, modZsize = 0.12, xw = 0.094, yh = 0.094, dist = 1.6)
+* 
+*  	MERLIN simulation; TS1 baseline model. 
+* 	In this example,  xw and yh are chosen to be identical to the shutter opening dimension.
+* 	dist = 1.6 is the real distance to the shutter front face:
+*       (This is TimeOffset value (=160 [cm]) from TS1verBase2016_LH8020_newVM-var_South04_Merlin.mcstas file.)
+* 
+*  N.B. Absolute normalization: The result of the Mc-Stas simulation will show neutron intensity for beam current of 1 uA.
+* 
+* 
+*  %P
+*  INPUT PARAMETERS:
+* 
+*  Face:     [string]  TS1 (or TS2) instrument McStas filename
+*  E0:          [meV] Lower edge of energy distribution
+*  E1:          [meV] Upper edge of energy distribution
+*  modXsize:      [m]   Moderator width
+*  modZsize:      [m]   Moderator height
+*  xw:            [m]   Width of focusing rectangle 
+*  yh:            [m]   Height of focusing rectangle
+*  dist:          [m]   Distance from moderator surface to the focusing rectangle
+*  verbose:     [int]     Flag to output debugging information
+*  beamcurrent: [uA]      ISIS beam current
+*  
+*  %E
+*******************************************************************************/
+DEFINE COMPONENT Commodus_I3
+SETTING PARAMETERS (string Face="TS1_S04_Merlin.mcstas",E0, E1, modPosition=0,
+  dist=1.7, int verbose=0, beamcurrent=1,
+  modXsize=0.12,modZsize=0.12,xw=0.094,yh=0.094)
+
+SHARE INHERIT ViewModISIS
+
+DECLARE INHERIT ViewModISIS EXTEND
+%{
+  double xwidth;
+  double yheight;
+  double focus_xw; 
+  double focus_yh;
+%}
+
+INITIALIZE
+%{
+  xwidth=modXsize;
+  yheight=modZsize;
+  focus_xw=xw;
+  focus_yh=yh;
+
+  /* READ IN THE ENERGY FILE */
+  FILE* Tfile;
+
+  Nsim=mcget_ncount();  // Number of points requested.
+
+  Tfile=openFile(Face);	              // Get open file
+  rtE0=convertEnergy(E0);
+  rtE1=convertEnergy(E1);
+  orderEnergy(&rtE0,&rtE1);
+
+  readHtable(Tfile,rtE0,rtE1, &TS, modPosition, xwidth, yheight, verbose);
+  fclose(Tfile);
+  // Below pragma was needed with PGI 19.x, compilation fails with NVC 20.7 
+  //#pragma acc declare create( TS )
+  /**********************************************************************/
+
+  Tnpts=0;
+  Ncount=0;
+
+  fprintf(stderr,"Face == %s \n",Face);
+  fprintf(stderr,"Number of Energy Points == %d\n",TS.nEnergy);
+  if (dist<0)
+    {
+      dist=TS.rdumMid;
+      fprintf(stderr,"Setting distance to moderator surface == %g\n",
+	      dist);
+    }
+  /* Do solid angle correction */
+  angleArea= strArea(TS, focus_xw, focus_yh, dist);
+  // Below pragma was needed with PGI 19.x, compilation fails with NVC 20.7 
+  //#pragma acc update host( TS )
+  fprintf(stderr,"Totals:: %g %d %d \n",TS.Total,TS.nEnergy,TS.nTime);
+
+
+%}
+
+TRACE INHERIT ViewModISIS
+
+MCDISPLAY INHERIT ViewModISIS
+
+END

mcstas-comps/contrib/ISIS_moderator.comp

@@ -127,6 +127,12 @@ double polInterp(double* X,double* Y,int Psize,double Aim)
   double out,errOut;         /* out put variables */
   double *C = malloc(Psize*sizeof(double));
   double *D = malloc(Psize*sizeof(double));
+  if (!C || !D) {
+    #ifndef OPENACC
+    fprintf(stderr,"Error in ISIS_moderator: memory allocation failure. Exit!\n");
+    exit(-1);
+    #endif
+  }
   double testDiff,diff;
 
   double w,den,ho,hp;           /* intermediate variables */
@@ -225,6 +231,10 @@ int cmdnumberD(char *mc,double* num)
 	(mc[i]=='\t' || mc[i]==' '  || mc[i]==',');i++);
   if(i==len || !mc[i]) return 0;
   ss=malloc(sizeof(char)*(len+1));
+  if(!ss) {
+    fprintf(stderr,"Error in ISIS_moderator: memory allocation failure. Exit!\n");
+    exit(-1);
+  }
 
   for(;i<len && mc[i]!='\n' && mc[i]
 	&& mc[i]!='\t' && mc[i]!=' ' && mc[i]!=',';i++)
@@ -392,10 +402,8 @@ int readHtable(FILE* TFile,const double Einit,const double Eend, Source *TS)
     Status Flag::
     Ftime=1 :: [time ] Reading Time : Data : Err [Exit on Total]
 
-
-    /*
     Double Read File to determine how many bins and
-    memery size
+    memory size
   */
   if (!TFile) return(0);
   Ea=0.0;

mcstas-comps/contrib/Multilayer_Sample.comp

@@ -10,7 +10,7 @@
 * Date: June 2010
 * Origin: ISIS
 *
-* Multilayer Reflecting sample using matrix Formula.
+* Multilayer Reflecting sample using matrix Formula, v2 with 'nrepeats' for repeating multilayer structure.
 *
 * %D
 *
@@ -45,12 +45,13 @@
 * target_index: [1] Used in combination with focus_xw and focus_yh to indicate solid angle for incoherent scattering.
 * focus_xw: [m] Used in combination with target_index and focus_yh to indicate solid angle for incoherent scattering.
 * focus_yh: [m] Used in combination with focus_xw and target_index to indicate solid angle for incoherent scattering.
+* nrepeats: [1] Number of repeats of the block of layers appart from the superphase and substrate
 *
 * %E
 *******************************************************************************/
 
 DEFINE COMPONENT Multilayer_Sample
-SETTING PARAMETERS (xwidth = 0.2, zlength = 0.2, nlayer=1,
+SETTING PARAMETERS (xwidth=0.2, zlength=0.2, int nlayer=1, int nrepeats=1,
   vector sldPar={0.0,2.0e-6,0.0e-6}, 
   vector dPar={20.0},
   vector sigmaPar={5.0,5.0}, 
@@ -174,6 +175,7 @@ TRACE
       double dbl0,dbl1,tvar;
       int arrsize=nlayer+2;
       int i;
+      int j;
 
       gsl_complex rnf,rnf1;
       gsl_complex a12t,a22t,cr,c0,ci;
@@ -252,7 +254,9 @@ TRACE
       gsl_matrix_complex_set(a1,1,1,cr);
 
       if(nlayer > 0){
-            for(i=1;i<nlayer+1;i++){
+        if(nrepeats > 1){
+          for(j=1;j<nrepeats-1;j++){
+            for(i=1;i<nlayer;i++){
 		btm=gsl_complex_mul(gsl_matrix_complex_get(betan,i,0),ci);
                 btm1=gsl_complex_mul(gsl_complex_mul_real(gsl_matrix_complex_get(betan,i,0),-1.0),ci);
 	        cbtm=gsl_complex_exp(btm);
@@ -284,7 +288,57 @@ TRACE
 		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,1));
                 gsl_matrix_complex_set(a3,0,1,gsl_complex_add(tcvar2,tcvar4));
 
-       		tcvar1=gsl_matrix_complex_get(a1,1,0);
+                tcvar1=gsl_matrix_complex_get(a1,1,0);
+                tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,0));
+                tcvar3=gsl_matrix_complex_get(a1,1,1);
+                tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,0));
+                gsl_matrix_complex_set(a3,1,0,gsl_complex_add(tcvar2,tcvar4));
+
+      		tcvar1=gsl_matrix_complex_get(a1,1,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,1));
+		tcvar3=gsl_matrix_complex_get(a1,1,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,1));
+                gsl_matrix_complex_set(a3,1,1,gsl_complex_add(tcvar2,tcvar4));
+
+		gsl_matrix_complex_set(a1,0,0,gsl_matrix_complex_get(a3,0,0));
+		gsl_matrix_complex_set(a1,0,1,gsl_matrix_complex_get(a3,0,1));
+		gsl_matrix_complex_set(a1,1,0,gsl_matrix_complex_get(a3,1,0));
+		gsl_matrix_complex_set(a1,1,1,gsl_matrix_complex_get(a3,1,1));
+            }
+          } // Exit the nrepeats block and do the final repeat with the substrate
+          for(i=1;i<nlayer+1;i++){
+		btm=gsl_complex_mul(gsl_matrix_complex_get(betan,i,0),ci);
+                btm1=gsl_complex_mul(gsl_complex_mul_real(gsl_matrix_complex_get(betan,i,0),-1.0),ci);
+		cbtm=gsl_complex_exp(btm);
+                cbtm1=gsl_complex_exp(btm1);
+    	        gsl_matrix_complex_set(a2,0,0,cbtm);
+        	tcvar1=gsl_matrix_complex_get(pfn,i,0);
+      		tcvar2=gsl_matrix_complex_get(pfn,i+1,0);
+      		if(GSL_REAL(gsl_complex_add(tcvar1,tcvar2))!=0.0 || GSL_IMAG(gsl_complex_add(tcvar1,tcvar2))!=0.0){
+		  a22t=gsl_complex_div(gsl_complex_sub(tcvar1,tcvar2),gsl_complex_add(tcvar1,tcvar2));
+      		} else {
+      		  a22t=c0;
+		}
+		tcvar3=gsl_complex_mul(tcvar1,tcvar2);
+      		tcvar4=gsl_complex_mul_real(tcvar3,-1.0*tlc*sigmaPar[i]*sigmaPar[i]);
+		a22t=gsl_complex_mul(a22t,gsl_complex_exp(tcvar4));
+      		gsl_matrix_complex_set(a2,0,1,gsl_complex_mul(a22t,cbtm));
+     		gsl_matrix_complex_set(a2,1,0,gsl_complex_mul(a22t,cbtm1));
+        	gsl_matrix_complex_set(a2,1,1,cbtm1);
+
+		tcvar1=gsl_matrix_complex_get(a1,0,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,0));
+		tcvar3=gsl_matrix_complex_get(a1,0,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,0));
+        	gsl_matrix_complex_set(a3,0,0,gsl_complex_add(tcvar2,tcvar4));
+
+		tcvar1=gsl_matrix_complex_get(a1,0,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,1));
+		tcvar3=gsl_matrix_complex_get(a1,0,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,1));
+        	gsl_matrix_complex_set(a3,0,1,gsl_complex_add(tcvar2,tcvar4));
+
+    		tcvar1=gsl_matrix_complex_get(a1,1,0);
 		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,0));
 		tcvar3=gsl_matrix_complex_get(a1,1,1);
 		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,0));
@@ -296,12 +350,63 @@ TRACE
 		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,1));
                 gsl_matrix_complex_set(a3,1,1,gsl_complex_add(tcvar2,tcvar4));
 
+		gsl_matrix_complex_set(a1,0,0,gsl_matrix_complex_get(a3,0,0));
+		gsl_matrix_complex_set(a1,0,1,gsl_matrix_complex_get(a3,0,1));
+		gsl_matrix_complex_set(a1,1,0,gsl_matrix_complex_get(a3,1,0));
+		gsl_matrix_complex_set(a1,1,1,gsl_matrix_complex_get(a3,1,1));
+          }
+	} else { // if we are not repeating then just do the same as the original
+          for(i=1;i<nlayer+1;i++){
+		btm=gsl_complex_mul(gsl_matrix_complex_get(betan,i,0),ci);
+                btm1=gsl_complex_mul(gsl_complex_mul_real(gsl_matrix_complex_get(betan,i,0),-1.0),ci);
+	        cbtm=gsl_complex_exp(btm);
+                cbtm1=gsl_complex_exp(btm1);
+    	        gsl_matrix_complex_set(a2,0,0,cbtm);
+        	tcvar1=gsl_matrix_complex_get(pfn,i,0);
+      		tcvar2=gsl_matrix_complex_get(pfn,i+1,0);
+      		if(GSL_REAL(gsl_complex_add(tcvar1,tcvar2))!=0.0 || GSL_IMAG(gsl_complex_add(tcvar1,tcvar2))!=0.0){
+		  a22t=gsl_complex_div(gsl_complex_sub(tcvar1,tcvar2),gsl_complex_add(tcvar1,tcvar2));
+      		} else {
+		  a22t=c0;
+		}
+		tcvar3=gsl_complex_mul(tcvar1,tcvar2);
+      		tcvar4=gsl_complex_mul_real(tcvar3,-1.0*tlc*sigmaPar[i]*sigmaPar[i]);
+		a22t=gsl_complex_mul(a22t,gsl_complex_exp(tcvar4));
+      		gsl_matrix_complex_set(a2,0,1,gsl_complex_mul(a22t,cbtm));
+     		gsl_matrix_complex_set(a2,1,0,gsl_complex_mul(a22t,cbtm1));
+        	gsl_matrix_complex_set(a2,1,1,cbtm1);
+
+		tcvar1=gsl_matrix_complex_get(a1,0,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,0));
+		tcvar3=gsl_matrix_complex_get(a1,0,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,0));
+       	        gsl_matrix_complex_set(a3,0,0,gsl_complex_add(tcvar2,tcvar4));
+
+		tcvar1=gsl_matrix_complex_get(a1,0,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,1));
+		tcvar3=gsl_matrix_complex_get(a1,0,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,1));
+       	        gsl_matrix_complex_set(a3,0,1,gsl_complex_add(tcvar2,tcvar4));
+
+	   	tcvar1=gsl_matrix_complex_get(a1,1,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,0));
+		tcvar3=gsl_matrix_complex_get(a1,1,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,0));
+    	        gsl_matrix_complex_set(a3,1,0,gsl_complex_add(tcvar2,tcvar4));
+
+		tcvar1=gsl_matrix_complex_get(a1,1,0);
+		tcvar2=gsl_complex_mul(tcvar1,gsl_matrix_complex_get(a2,0,1));
+		tcvar3=gsl_matrix_complex_get(a1,1,1);
+		tcvar4=gsl_complex_mul(tcvar3,gsl_matrix_complex_get(a2,1,1));
+    	        gsl_matrix_complex_set(a3,1,1,gsl_complex_add(tcvar2,tcvar4));
+
 		gsl_matrix_complex_set(a1,0,0,gsl_matrix_complex_get(a3,0,0));
 		gsl_matrix_complex_set(a1,0,1,gsl_matrix_complex_get(a3,0,1));
 		gsl_matrix_complex_set(a1,1,0,gsl_matrix_complex_get(a3,1,0));
 		gsl_matrix_complex_set(a1,1,1,gsl_matrix_complex_get(a3,1,1));
             }
-      }
+         }
+      }        
       ac1=gsl_matrix_complex_get(a1,1,0);
       ac2=gsl_complex_conjugate(ac1);
       ac3=gsl_matrix_complex_get(a1,0,0);