kw_monitor
Farhi
farhi at ill.fr
Wed Oct 6 20:14:11 CEST 1999
Hy Kristian
Here are the files for the Optimizer.
Yes I use the in14_6 instrument to evaluate the optimizer efficiency.
You can plot the kw resolution function at the end.
as a reference I also include a in14_6n instrument without optimization.
If you got some pbs in compiling, tell me...
here follows a Log for a simulation :
mcstas in14_6.instr
gcc -O2 -o in14_6 in14_6.c -lm
farhi at macfarhi Instr:153> in14_6 -n 1e6 KI=2.66 WN=0.03 ORDER=1 MHV=3
Instrument : IN14, v60 (21/09/99) on macfarhi.
Flat source, m=3.00 noze, width 0.03
Monochromator : (DM = 3.355)
A1 = 20.61, A2 = 41.22 (deg)
Ki = 2.66 Angs-1 (Energy = 14.7 meV, Velocity = 1.68e+03 m/s)
RM = 0.669 Deg, RH = 0.714 Deg
Sample :
A3 = 0.00, A4 = 0.00 (deg)
Energy transfert 0 meV, Moment transfert 0 Angs-1
Analyser : (DA = 3.355)
A5 = 20.61, A6 = 41.22 (deg)
Kf = 2.66 Angs-1 (Energy = 14.7 meV, Velocity = 1.68e+03 m/s)
RA = 1.21 Deg
Detectors :
>> OPTIM_PHASE_SET_REF (1000 neutrons)
>> AUTO monitor has reached 100 counts (non optimized, absorbed)
0.119690 119690
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons) from Ref
Counts : reference = 119690, passing = 119690, monitor = 100
Flux : reference = 1.9e+11, passing = 1.9e+11, monitor = 3.7e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 106379
Counts : reference = 119690, passing = 119690, monitor = 4539
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.4e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 212824
Counts : reference = 119690, passing = 119690, monitor = 4544
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.4e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 319299
Counts : reference = 119690, passing = 119690, monitor = 4580
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.4e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 425880
Counts : reference = 119690, passing = 119690, monitor = 4588
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.1e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 532178
Counts : reference = 119690, passing = 119690, monitor = 4606
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.4e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 638681
Counts : reference = 119690, passing = 119690, monitor = 4629
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.3e+07
>> OPTIM_PHASE_SET_SOURCE (119690 neutrons)
Number of redirections : 745164
Counts : reference = 119690, passing = 119690, monitor = 4658
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 3.6e+07
End of optimization
Optim_Normal_Monitor_Counts = 100 (2 steps), Optim_Total_Monitor_Counts =
32968
Optimizer speedup : 78.9
Number of redirections : 761697
Counts : reference = 119690, passing = 28691, monitor = 724
Flux : reference = 1.9e+11, passing = 1.9e+12, monitor = 6.5e+06
Detector: kw1_I=7.07511e+11 kw1_ERR=1.34359e+10 kw1_N=60746
Detector: PSDSample_I=2.41983e+08 PSDSample_ERR=2.85056e+06 PSDSample_N=32868
Detector: Int4Sample_I=2.3616e+08 Int4Sample_ERR=2.81932e+06
Int4Sample_N=29140
Detector: Flux1Sample_I=3.57466e+07 Flux1Sample_ERR=1.05143e+06
Flux1Sample_N=3677
Detector: ESample_I=3.54956e+07 ESample_ERR=1.04524e+06 ESample_N=3654
Detector: DivSample_I=2.9073e+07 DivSample_ERR=950905 DivSample_N=2399
Detector: kw_I=3.56749e+07 kw_ERR=1.04785e+06 kw_N=3669
Output files : sim/i60_k27w30d1m3.psd sim/i60_k27w30d1m3.nrj
sim/i60_k27w30d1m3.div sim/i60_k27w30d1m3.kw
Cheers !
--
Emmanuel FARHI, http://www.ill.fr/tas/people/Farhi.html \|/ ____ \|/
TAS-Group, Institut Laue-Langevin (ILL) Grenoble ~@-/ oO \-@~
Avenue des Martyrs, BP 156, 38042 Grenoble Cedex 9,France /_( \__/ )_\
Work :Tel (33/0) 4 76 20 71 83. Fax (33/0) 4 76 48 39 06 \__U_/
La Grande Arche, Chateau d'Uriage, 38410 Saint Martin d'Uriage 04 76 59 73 94
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/19991006/7c9b9cb2/attachment.html>
-------------- next part --------------
/*******************************************************************************
*
* McStas, version 1.1, released
* Maintained by Kristian Nielsen and Kim Lefmann,
* Risoe National Laboratory, Roskilde, Denmark
*
* Component: Source_Optimizer
*
* Written by: EF, 17 Sept 1999
*
* Usage: A component that optimizes the neutron flux passing through the
* Source_optimizer in order to have the maximum flux at the
* Monitor_Optimizer position.
*
* Principle: The optimizer first computes neutron state parameter limits
* (step 1), and then records a Reference source (step 2). The optimization then
* starts (step 3), and focuses new neutrons on the Monitor_Optimizer.
*
* Options: The optimized source can be computed regularly ('continuous' option)
* or only once ('fixed'). The energy distribution can be kept during optimization
* ('keepE') or released ('freeE'). The time spent in steps 1 and 2 can be reduced
* for a better optimization ('auto'). The neutrons passing during steps 1 and 2
* can be absorbed for a better neutron weight distribution.
*
* Source_optimizer should be placed just after the source.
* Monitor_Optimizer should be placed at position to optimize.
*
* INPUT PARAMETERS:
*
* xmin: Lower x bound of optimizer opening (m)
* xmax: Upper x bound of optimizer opening (m)
* ymin: Lower y bound of optimizer opening (m)
* ymax: Upper y bound of optimizer opening (m)
* bins: Number of cells for sampling of neutron state distribution (10 is default)
* step: Optimizer step (%, 10 is default)
* The two first steps are not optimized.
* keep: Percentage of initial source distribution kept (%, 10 is default)
* file: Filename where to save optimized source distributions
* options: string that can contain
* 'continuous' for continuous source optimization (default)
* 'fixed' same as 'not continuous' optimization
* 'keepE' to keep energy and velocity if pos. (default)
* 'freeE' same as 'no keepE'
* 'verbose' displays optimization process
* 'auto' uses the shortest possible 'step 1' and 'step 2'
* and set 'step' as required. Make sure to have
* enough neutrons in simulation (ncount)
* 'absorb' absorbs step 1 and 2 neutrons if possible.
*
* parameters bins, step, keep can be -1 for default values.
*
* OUTPUT PARAMETERS:
*
* distributions if filename is not empty ("")
*
* EXAMPLE: I use the following settings
* xmin = -0.03, xmax = 0.03,
* ymin = -0.06, ymax = 0.06,
* bins = -1, step = -1, keep = -1,
* file="source.optim",
* options="absorb+auto"
*
*******************************************************************************/
/* History :
Sep 17 1999 : v0.00 first release (not effective)
Sep 26 1999 : v0.01 New_Source for continuous optimisation
Sep 27 1999 : optimizer is ok, but not efficient
Sep 29 1999 : v0.02 re-direct 'bad' neutrons instead of ABSORB (rand generator for nothing)
*/
DEFINE COMPONENT Source_Optimizer
DEFINITION PARAMETERS (xmin, xmax, ymin, ymax, bins, step, keep,file,options)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
#ifndef MCSTAS_GEN_OPTIM /* McStas General optimizer ID */
#define MCSTAS_GEN_OPTIM
#else
#error McStas : Source_Optimizer should only be used once per instrument
#endif
#ifndef FLT_MAX
#define FLT_MAX 1e37
#endif
#define OPTIM_PHASE_SET_LIMITS 0 /* set array limits to 0, then ask for GET_LIMITS */
#define OPTIM_PHASE_GET_LIMITS 1 /* compute array limits, then ask for SET_REF */
#define OPTIM_PHASE_SET_REF 2 /* set Ref and New_Source to to 0, then ask for GET_REF */
#define OPTIM_PHASE_GET_REF 3 /* compute Ref (and New_Source in Monitor), then ask for SET_SOURCE */
#define OPTIM_PHASE_SET_SOURCE 4 /* set Source to Ref*x%+New_Source, normalize to Ref, Passing to 0, then ask for OPTIM */
#define OPTIM_PHASE_OPTIM 5 /* Optimize and get New_Source (continuous optimization), then reask SET_SOURCE when required */
#define OPTIM_MOD_X 1
#define OPTIM_MOD_Y 2
#define OPTIM_MOD_VX 4
#define OPTIM_MOD_VY 8
#define OPTIM_MOD_VZ 16
#define OPTIM_MOD_S1 32
#define OPTIM_MOD_S2 64
#define OPTIM_MOD_P 128
/* These are distribution arrays[bins] within limits
* flux is kept during optimisation
* NOT stored : z is the position of this component
* t time (linked to z)
*/
/* initial Reference distribution arrays (for weights) */
long *Optim_Reference_x;
long *Optim_Reference_y;
long *Optim_Reference_vx;
long *Optim_Reference_vy;
long *Optim_Reference_vz;
long *Optim_Reference_s1;
long *Optim_Reference_s2;
long *Optim_Reference_p;
/* optimized Source distribution arrays (to reach) */
long *Optim_Source_x;
long *Optim_Source_y;
long *Optim_Source_vx;
long *Optim_Source_vy;
long *Optim_Source_vz;
long *Optim_Source_s1;
long *Optim_Source_s2;
long *Optim_Source_p;
/* optimized New_Source distribution arrays (to reach in next step, passed to Source) */
long *Optim_New_Source_x;
long *Optim_New_Source_y;
long *Optim_New_Source_vx;
long *Optim_New_Source_vy;
long *Optim_New_Source_vz;
long *Optim_New_Source_s1;
long *Optim_New_Source_s2;
long *Optim_New_Source_p;
/* Passing distribution arrays (should grow to reach Source) */
long *Optim_Passing_x;
long *Optim_Passing_y;
long *Optim_Passing_vx;
long *Optim_Passing_vy;
long *Optim_Passing_vz;
long *Optim_Passing_s1;
long *Optim_Passing_s2;
long *Optim_Passing_p;
/* limits for state parameters */
/* x and y are Optimizer dimensions (input parameters) */
double Optim_x_min, Optim_x_max;
double Optim_y_min, Optim_y_max;
double Optim_vx_min, Optim_vx_max;
double Optim_vy_min, Optim_vy_max;
double Optim_vz_min, Optim_vz_max;
double Optim_s1_min, Optim_s1_max;
double Optim_s2_min, Optim_s2_max;
double Optim_p_min, Optim_p_max;
int Optim_index=0; /* a running Optim_index */
int Optim_index_x=0; /* indexes for last neutron that passed through */
int Optim_index_y=0;
int Optim_index_vx=0;
int Optim_index_vy=0;
int Optim_index_vz=0;
int Optim_index_s1=0;
int Optim_index_s2=0;
int Optim_index_p=0;
int Optim_good_x=0; /* indexes for last 'good' neutron that passed through */
int Optim_good_y=0;
int Optim_good_vx=0;
int Optim_good_vy=0;
int Optim_good_vz=0;
int Optim_good_s1=0;
int Optim_good_s2=0;
int Optim_good_p=0;
int Optim_bins;
long Optim_n_redirect; /* number of consecutive ABSORB */
int Optim_Phase; /* Optimizer function */
long Optim_Phase_Counts =0; /* neutron counts to achieve in each Phase */
long Optim_Phase_Counts_L =0; /* neutron counts to achieve in Limits Phase */
long Optim_Phase_Counts_R =0; /* neutron counts to achieve in Reference Phase */
char Optim_Flag_Continuous =0; /* 1 : continuous Source optimization */
char Optim_Flag_Recycle =0; /* record of neutron state changes by OPTIM_MOD_xx */
char Optim_Flag_KeepE =0; /* 1 : keep E as poss. when recycling */
/* i.e. keep E and |v| distribution */
char Optim_Flag_Verbose =0; /* displays optimization informations */
char Optim_Flag_Absorb =0; /* 1 means that first steps non optimized neutrons are absorbed */
char Optim_Flag_Auto =0; /* 1 is for minimum counts in 2 first steps */
long Optim_Limits_Counts =0; /* passing neutron counts in each Optim_Phase */
long Optim_Reference_Counts =0;
long Optim_Passing_Counts =0;
long Optim_Monitor_Counts =0;
double Optim_Limits_Flux =0; /* passing neutron flux in each Optim_Phase */
double Optim_Reference_Flux=0;
double Optim_Passing_Flux =0;
double Optim_Monitor_Flux =0;
double Optim_Monitor_pmax =0;
float Optim_keep;
float Optim_step;
double Optim_vx; /* save neutron characteristics for Monitor and ABSORDed->Redirected neutrons */
double Optim_vy;
double Optim_vz;
double Optim_x;
double Optim_y;
double Optim_s1;
double Optim_s2;
double Optim_p;
double Optim_v2;
double Optim_t1; /* tempo vars */
double Optim_t2;
double Optim_t3;
double Optim_u1; /* tempo vars */
double Optim_u2;
double Optim_u3;
int Optim_i1; /* tempo vars */
int Optim_i2;
int Optim_i3;
long Optim_Normal_Monitor_Counts = 0; /* counts without optim */
long Optim_Total_Monitor_Counts =0; /* final monitor counts */
FILE *hfile;
/* end declare */
%}
INITIALIZE
%{
Optim_n_redirect = 0;
Optim_Phase = OPTIM_PHASE_SET_LIMITS;
Optim_x_min = xmin;
Optim_x_max = xmax;
Optim_y_min = ymin;
Optim_y_max = ymax;
Optim_bins = (int)bins;
Optim_step = step;
Optim_keep = keep;
if (Optim_step < 0) Optim_step = .1; /* default values if -1 is given */
if (Optim_bins < 0) Optim_bins = 10;
if (Optim_keep < 0) Optim_keep = .1;
if (Optim_step >= 1) Optim_step = Optim_step/100; /* in case user gives % in 1-100 */
if (Optim_step < .01) Optim_step = .01;
if (Optim_step > 1) Optim_step = 1;
if (Optim_keep >= 1) Optim_keep = Optim_keep/100; /* in case user gives % in 1-100 */
if (Optim_keep < .01) Optim_keep = .01;
if (Optim_keep > 1) Optim_keep = 1;
Optim_Phase_Counts = mcget_ncount() * Optim_step;
Optim_Phase_Counts_L = Optim_Phase_Counts;
Optim_Phase_Counts_R = Optim_Phase_Counts;
if (Optim_bins < 1) Optim_bins = 1;
if (Optim_bins > 100) Optim_bins = 100;
if (strstr(options,"continuous")) Optim_Flag_Continuous = 1;
if (strstr(options,"fixed")) Optim_Flag_Continuous = 0;
if (strstr(options,"keepE")) Optim_Flag_KeepE = 1;
if (strstr(options,"freeE")) Optim_Flag_KeepE = 0;
if (strstr(options,"verbose")) Optim_Flag_Verbose = 1;
if (strstr(options,"auto"))
{
Optim_Flag_Auto = 1;
if (Optim_bins*10 < Optim_Phase_Counts) Optim_Phase_Counts_L = Optim_bins*100; /* need at least 10 counts per bin for Limits */
Optim_Phase_Counts_R = mcget_ncount();
Optim_Phase_Counts = mcget_ncount();
}
if (strstr(options,"absorb")) Optim_Flag_Absorb = 1;
if ((Optim_Source_x = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_y = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_vx = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_vy = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_vz = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_s1 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_s2 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Source_p = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_x = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_y = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_vx = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_vy = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_vz = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_s1 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_s2 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_New_Source_p = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_x = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_y = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_vx = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_vy = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_vz = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_s1 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_s2 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Passing_p = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_x = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_y = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_vx = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_vy = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_vz = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_s1 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_s2 = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
if ((Optim_Reference_p = (long*)malloc(Optim_bins * sizeof(long*))) == NULL) { fprintf(stderr,"Optimizer : not enough memory\n"); exit(-1); }
/* end initialize */
%}
TRACE
%{
PROP_Z0;
if (x>Optim_x_min && x<Optim_x_max && y>Optim_y_min && y<Optim_y_max)
{
Optim_vx = vx; /* save neutron characteristics for Monitor */
Optim_vy = vy;
Optim_vz = vz;
Optim_x = x;
Optim_y = y;
Optim_s1 = s1;
Optim_s2 = s2;
Optim_p = p;
Optim_v2 = (vx*vx+vy*vy+vz*vz); /* squared velocity */
/* handle Phase sequence */
if ((Optim_Phase == OPTIM_PHASE_GET_LIMITS)
&& (Optim_Limits_Counts >= Optim_Phase_Counts_L))
{
Optim_Phase = OPTIM_PHASE_SET_REF;
if (Optim_Flag_Verbose) printf(">> OPTIM_PHASE_SET_REF (%i neutrons)\n", Optim_Limits_Counts);
}
if ((Optim_Phase == OPTIM_PHASE_GET_REF)
&& (Optim_Reference_Counts >= Optim_Phase_Counts_R))
{
Optim_Phase = OPTIM_PHASE_SET_SOURCE;
Optim_Phase_Counts_R = Optim_Phase_Counts;
if (Optim_Flag_Verbose)
{
printf(">> OPTIM_PHASE_SET_SOURCE (%i neutrons) from Ref\n", Optim_Reference_Counts);
printf("Counts : reference = %i, passing = %i, monitor = %i\n", Optim_Reference_Counts, Optim_Reference_Counts, Optim_Monitor_Counts);
printf("Flux : reference = %.2g, passing = %.2g, monitor = %.2g\n", Optim_Reference_Flux, Optim_Reference_Flux, Optim_Monitor_Flux);
}
}
if ((Optim_Phase == OPTIM_PHASE_OPTIM)
&& (Optim_Passing_Counts >= Optim_Phase_Counts))
{
Optim_Phase = OPTIM_PHASE_SET_SOURCE;
if (Optim_Flag_Verbose)
{
printf(">> OPTIM_PHASE_SET_SOURCE (%i neutrons)\n", Optim_Passing_Counts);
printf("Number of redirections : %i\n",Optim_n_redirect);
printf("Counts : reference = %i, passing = %i, monitor = %i\n", Optim_Reference_Counts, Optim_Passing_Counts, Optim_Monitor_Counts);
printf("Flux : reference = %.2g, passing = %.2g, monitor = %.2g\n", Optim_Reference_Flux, Optim_Passing_Flux, Optim_Monitor_Flux);
}
}
/* handle Optim_Phase functions */
if (Optim_Phase == OPTIM_PHASE_SET_LIMITS) /* init : need to compute limits and flux */
{
Optim_Limits_Counts = 0;
Optim_Limits_Flux = 0;
Optim_vx_min = FLT_MAX; Optim_vx_max = -FLT_MAX;
Optim_vy_min = FLT_MAX; Optim_vy_max = -FLT_MAX;
Optim_vz_min = FLT_MAX; Optim_vz_max = -FLT_MAX;
Optim_s1_min = FLT_MAX; Optim_s1_max = -FLT_MAX;
Optim_s2_min = FLT_MAX; Optim_s2_max = -FLT_MAX;
Optim_p_min = FLT_MAX; Optim_p_max = -FLT_MAX;
Optim_Phase = OPTIM_PHASE_GET_LIMITS;
} /* end OPTIM_PHASE_SET_LIMITS */
if (Optim_Phase == OPTIM_PHASE_GET_LIMITS) /* init : need to compute limits and flux */
{
Optim_Limits_Counts++;
Optim_Limits_Flux += p;
if (x < Optim_x_min) Optim_x_min = x;
if (y < Optim_y_min) Optim_y_min = y;
if (x > Optim_x_max) Optim_x_max = x;
if (y > Optim_y_max) Optim_y_max = y;
if (vx < Optim_vx_min) Optim_vx_min = vx;
if (vx > Optim_vx_max) Optim_vx_max = vx;
if (vy < Optim_vy_min) Optim_vy_min = vy;
if (vy > Optim_vy_max) Optim_vy_max = vy;
if (vz < Optim_vz_min) Optim_vz_min = vz;
if (vz > Optim_vz_max) Optim_vz_max = vz;
if (p < Optim_p_min) Optim_p_min = p;
if (p > Optim_p_max) Optim_p_max = p;
if (s1 < Optim_s1_min) Optim_s1_min = s1;
if (s1 > Optim_s1_max) Optim_s1_max = s1;
if (s2 < Optim_s2_min) Optim_s2_min = s2;
if (s2 > Optim_s2_max) Optim_s2_max = s2;
if (Optim_Flag_Absorb) ABSORB;
} /* end if OPTIM_PHASE_GET_LIMITS */
if (Optim_Phase == OPTIM_PHASE_SET_REF) /* Set Ref and New_Source to 0 */
{
Optim_Reference_Counts = 0;
Optim_Reference_Flux = 0;
Optim_Monitor_Counts = 0; /* also counted as New_Source */
Optim_Monitor_Flux = 0;
for (Optim_index=0; Optim_index < Optim_bins; Optim_index++)
{
Optim_Reference_x[Optim_index] = 0; /* initial distribution will be recorded first */
Optim_Reference_y[Optim_index] = 0;
Optim_Reference_vx[Optim_index] = 0;
Optim_Reference_vy[Optim_index] = 0;
Optim_Reference_vz[Optim_index] = 0;
Optim_Reference_s1[Optim_index] = 0;
Optim_Reference_s2[Optim_index] = 0;
Optim_Reference_p[Optim_index] = 0;
Optim_New_Source_x[Optim_index] = 0; /* Monitor_Optimizer will compute the */
Optim_New_Source_y[Optim_index] = 0; /* optimized New_Source distribution */
Optim_New_Source_vx[Optim_index] = 0; /* that will become Source for Optim Optim_step */
Optim_New_Source_vy[Optim_index] = 0;
Optim_New_Source_vz[Optim_index] = 0;
Optim_New_Source_s1[Optim_index] = 0;
Optim_New_Source_s2[Optim_index] = 0;
Optim_New_Source_p[Optim_index] = 0;
} /* end for */
Optim_Phase = OPTIM_PHASE_GET_REF;
} /* end OPTIM_PHASE_SET_REF */
if (Optim_Phase == OPTIM_PHASE_GET_REF) /* now build the Reference in limits */
{ /* New_Source is set by Monitor_Optimizer */
Optim_Reference_Counts++;
Optim_Reference_Flux += p;
if (Optim_vx_max-Optim_vx_min)
Optim_index = (int)rint(Optim_bins * (vx -Optim_vx_min)/(Optim_vx_max-Optim_vx_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_vx[Optim_index]++;
if (Optim_vy_max-Optim_vy_min)
Optim_index = (int)rint(Optim_bins * (vy -Optim_vy_min)/(Optim_vy_max-Optim_vy_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_vy[Optim_index]++;
if (Optim_vz_max-Optim_vz_min)
Optim_index = (int)rint(Optim_bins * (vz -Optim_vz_min)/(Optim_vz_max-Optim_vz_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_vz[Optim_index]++;
if (Optim_x_max-Optim_x_min)
Optim_index = (int)rint(Optim_bins * (x -Optim_x_min)/(Optim_x_max-Optim_x_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_x[Optim_index]++;
if (Optim_y_max-Optim_y_min)
Optim_index = (int)rint(Optim_bins * (y -Optim_y_min)/(Optim_y_max-Optim_y_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_y[Optim_index]++;
if (Optim_s1_max-Optim_s1_min)
Optim_index = (int)rint(Optim_bins * (s1 -Optim_s1_min)/(Optim_s1_max-Optim_s1_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_s1[Optim_index]++;
if (Optim_s2_max-Optim_s2_min)
Optim_index = (int)rint(Optim_bins * (s2 -Optim_s2_min)/(Optim_s2_max-Optim_s2_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_s2[Optim_index]++;
if (Optim_p_max-Optim_p_min)
Optim_index = (int)rint(Optim_bins * (p -Optim_p_min)/(Optim_p_max-Optim_p_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_Reference_p[Optim_index]++;
} /* end if OPTIM_PHASE_GET_REF */
if (Optim_Phase == OPTIM_PHASE_SET_SOURCE) /* Define optimized Source (normalized to Reference) */
{
if (Optim_Monitor_Counts)
Optim_t1 = (1 - Optim_keep) * Optim_Reference_Counts/Optim_Monitor_Counts;
else
Optim_t1 = 0;
Optim_Passing_Counts = 0;
Optim_Passing_Flux = 0;
if (Optim_Normal_Monitor_Counts == 0) Optim_Normal_Monitor_Counts = Optim_Total_Monitor_Counts; /* 2 first un-optimized steps */
Optim_Monitor_Counts = 0; /* also counted as New_Source */
Optim_Monitor_Flux = 0;
for (Optim_index = 0; Optim_index < Optim_bins; Optim_index++)
{ /* get Optim_keep % of Reference, and 1-Optim_keep% of New_Source normalized to Reference Counts */
if (Optim_Flag_Continuous | (Optim_n_redirect == 0))
{
Optim_Source_x[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_x[Optim_index]) + Optim_t1 * Optim_New_Source_x[Optim_index] );
Optim_Source_y[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_y[Optim_index]) + Optim_t1 * Optim_New_Source_y[Optim_index] );
Optim_Source_vx[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_vx[Optim_index]) + Optim_t1 * Optim_New_Source_vx[Optim_index] );
Optim_Source_vy[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_vy[Optim_index]) + Optim_t1 * Optim_New_Source_vy[Optim_index] );
Optim_Source_vz[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_vz[Optim_index]) + Optim_t1 * Optim_New_Source_vz[Optim_index] );
Optim_Source_s1[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_s1[Optim_index]) + Optim_t1 * Optim_New_Source_s1[Optim_index] );
Optim_Source_s2[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_s2[Optim_index]) + Optim_t1 * Optim_New_Source_s2[Optim_index] );
Optim_Source_p[Optim_index] = (long)(rint(Optim_keep * Optim_Reference_p[Optim_index]) + Optim_t1 * Optim_New_Source_p[Optim_index] );
if (Optim_New_Source_x[Optim_index] > Optim_New_Source_x[Optim_good_x]) Optim_good_x = Optim_index;
if (Optim_New_Source_y[Optim_index] > Optim_New_Source_y[Optim_good_y]) Optim_good_y = Optim_index;
if (Optim_New_Source_vx[Optim_index] > Optim_New_Source_vx[Optim_good_vx]) Optim_good_vx = Optim_index;
if (Optim_New_Source_vy[Optim_index] > Optim_New_Source_vy[Optim_good_vy]) Optim_good_vy = Optim_index;
if (Optim_New_Source_vz[Optim_index] > Optim_New_Source_vz[Optim_good_vz]) Optim_good_vz = Optim_index;
if (Optim_New_Source_s1[Optim_index] > Optim_New_Source_s1[Optim_good_s1]) Optim_good_s1 = Optim_index;
if (Optim_New_Source_s2[Optim_index] > Optim_New_Source_s2[Optim_good_s2]) Optim_good_s2 = Optim_index;
if (Optim_New_Source_p[Optim_index] > Optim_New_Source_p[Optim_good_p]) Optim_good_p = Optim_index;
}
Optim_Passing_x[Optim_index] = 0; /* Passing neutrons will then reach Source */
Optim_Passing_y[Optim_index] = 0; /* weights will be adapted to match Reference */
Optim_Passing_vx[Optim_index] = 0;
Optim_Passing_vy[Optim_index] = 0;
Optim_Passing_vz[Optim_index] = 0;
Optim_Passing_s1[Optim_index] = 0;
Optim_Passing_s2[Optim_index] = 0;
Optim_Passing_p[Optim_index] = 0;
Optim_New_Source_x[Optim_index] = 0; /* Init of next Source */
Optim_New_Source_y[Optim_index] = 0;
Optim_New_Source_vx[Optim_index] = 0;
Optim_New_Source_vy[Optim_index] = 0;
Optim_New_Source_vz[Optim_index] = 0;
Optim_New_Source_s1[Optim_index] = 0;
Optim_New_Source_s2[Optim_index] = 0;
Optim_New_Source_p[Optim_index] = 0;
} /* end for */
Optim_Phase = OPTIM_PHASE_OPTIM;
} /* end OPTIM_PHASE_SET_SOURCE */
if (Optim_Phase == OPTIM_PHASE_OPTIM) /* Use optimized Source */
{
Optim_Flag_Recycle = 0;
Optim_index_x = Optim_good_x;
Optim_index_y = Optim_good_y;
Optim_index_vx= Optim_good_vx;
Optim_index_vy= Optim_good_vy;
Optim_index_vz= Optim_good_vz;
Optim_index_s1= Optim_good_s1;
Optim_index_s2= Optim_good_s2;
Optim_index_p = Optim_good_p;
if (Optim_vz_max-Optim_vz_min)
Optim_index = (int)rint(Optim_bins * (vz -Optim_vz_min)/(Optim_vz_max-Optim_vz_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_vz[Optim_index] >= Optim_Source_vz[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_VX;
Optim_vz += (Optim_index_vz-Optim_index)*(Optim_vz_max - Optim_vz_min)/Optim_bins;
}
else
Optim_index_vz = Optim_index;
if (Optim_vx_max-Optim_vx_min)
Optim_index = (int)rint(Optim_bins * (vx -Optim_vx_min)/(Optim_vx_max-Optim_vx_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_vx[Optim_index] >= Optim_Source_vx[Optim_index]) /* distribution achieved : redirect neutron near last neutron characteristic */
{
Optim_Flag_Recycle |= OPTIM_MOD_VY;
Optim_vx += (Optim_index_vx-Optim_index)*(Optim_vx_max - Optim_vx_min)/Optim_bins;
}
else
Optim_index_vx = Optim_index;
if (Optim_vy_max-Optim_vy_min)
Optim_index = (int)rint(Optim_bins * (vy -Optim_vy_min)/(Optim_vy_max-Optim_vy_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_vy[Optim_index] >= Optim_Source_vy[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_VZ;
Optim_vy += (Optim_index_vy-Optim_index)*(Optim_vy_max - Optim_vy_min)/Optim_bins;
}
else
Optim_index_vy = Optim_index;
if ((Optim_Flag_Recycle & (OPTIM_MOD_VX|OPTIM_MOD_VY|OPTIM_MOD_VZ)) && Optim_Flag_KeepE)
{ /* now try to keep E distribution */
Optim_t1 = Optim_v2 - Optim_vz*Optim_vz - Optim_vy*Optim_vy;
Optim_t2 = Optim_v2 - Optim_vz*Optim_vz - Optim_vx*Optim_vx;
Optim_t3 = Optim_v2 - Optim_vx*Optim_vx - Optim_vy*Optim_vy;
/* we affect the component wich is the most optimized (largest Source/Ref) */
if ((Optim_vx_max-Optim_vx_min) && (Optim_t1 > 0))
{
Optim_t1 = sqrt(Optim_t1);
if (vx < 0) Optim_t1 = -Optim_t1;
Optim_i1 = (int)rint(Optim_bins * (Optim_t1 -Optim_vx_min)/(Optim_vx_max-Optim_vx_min));
if (Optim_i1 < 0) Optim_i1 = 0;
if (Optim_i1 >= Optim_bins) Optim_i1 = Optim_bins - 1;
Optim_u1 = (double)Optim_Source_vx[Optim_i1]/(Optim_Reference_vx[Optim_i1]+1);
}
else
Optim_u1 = 0;
if ((Optim_vy_max-Optim_vy_min) && (Optim_t2 > 0))
{
Optim_t2 = sqrt(Optim_t2);
if (vy < 0) Optim_t2 = -Optim_t2;
Optim_i2 = (int)rint(Optim_bins * (Optim_t2 -Optim_vy_min)/(Optim_vy_max-Optim_vy_min));
if (Optim_i2 < 0) Optim_i2 = 0;
if (Optim_i2 >= Optim_bins) Optim_i2 = Optim_bins - 1;
Optim_u2 = (double)Optim_Source_vy[Optim_i2]/(Optim_Reference_vy[Optim_i2]+1);
}
else
Optim_u2 = 0;
if ((Optim_vz_max-Optim_vz_min) && (Optim_t3 > 0))
{
Optim_t3 = sqrt(Optim_t3);
if (vz < 0) Optim_t3 = -Optim_t3;
Optim_i3 = (int)rint(Optim_bins * (Optim_t3 -Optim_vz_min)/(Optim_vz_max-Optim_vz_min));
if (Optim_i3 < 0) Optim_i3 = 0;
if (Optim_i3 >= Optim_bins) Optim_i3 = Optim_bins - 1;
Optim_u3 = (double)Optim_Source_vz[Optim_i3]/(Optim_Reference_vz[Optim_i3]+1);
}
else
Optim_u3 = 0;
if ((Optim_u1 > Optim_u2) && (Optim_u1 > Optim_u3))
{
Optim_vx = Optim_t1;
Optim_index_vx = Optim_i1;
Optim_Flag_Recycle |= OPTIM_MOD_VX;
Optim_index = -1;
}
if ((Optim_u2 > Optim_u1) && (Optim_u2 > Optim_u3) )
{
Optim_vy = Optim_t2;
Optim_index_vy = Optim_i2;
Optim_Flag_Recycle |= OPTIM_MOD_VY;
Optim_index = -1;
}
if ((Optim_u3 > Optim_u1) && (Optim_u3 > Optim_u1))
{
Optim_vz = Optim_t3;
Optim_index_vz = Optim_i3;
Optim_Flag_Recycle |= OPTIM_MOD_VZ;
Optim_index = -1;
}
}
vx = Optim_vx;
vy = Optim_vy;
vz = Optim_vz;
if (Optim_x_max-Optim_x_min)
Optim_index = (int)rint(Optim_bins * (x -Optim_x_min)/(Optim_x_max-Optim_x_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_x[Optim_index] >= Optim_Source_x[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_X;
Optim_x += (Optim_index_x-Optim_index)*(Optim_x_max - Optim_x_min)/Optim_bins;
x = Optim_x;
}
else
Optim_index_x = Optim_index;
if (Optim_y_max-Optim_y_min)
Optim_index = (int)rint(Optim_bins * (y -Optim_y_min)/(Optim_y_max-Optim_y_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_y[Optim_index] >= Optim_Source_y[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_Y;
Optim_y += (Optim_index_y-Optim_index)*(Optim_y_max - Optim_y_min)/Optim_bins;
y = Optim_y;
}
else
Optim_index_y = Optim_index;
if (Optim_s1_max-Optim_s1_min)
Optim_index = (int)rint(Optim_bins * (s1 -Optim_s1_min)/(Optim_s1_max-Optim_s1_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_s1[Optim_index] >= Optim_Source_s1[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_S1;
Optim_s1 += (Optim_index_s1-Optim_index)*(Optim_s1_max - Optim_s1_min)/Optim_bins;
s1 = Optim_s1;
}
else
Optim_index_s1 = Optim_index;
if (Optim_s2_max-Optim_s2_min)
Optim_index = (int)rint(Optim_bins * (s2 -Optim_s2_min)/(Optim_s2_max-Optim_s2_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_s2[Optim_index] >= Optim_Source_s2[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_S2;
Optim_s2 += (Optim_index_s2-Optim_index)*(Optim_s2_max - Optim_s2_min)/Optim_bins;
s2 = Optim_s2;
}
else
Optim_index_s2 = Optim_index;
if (Optim_p_max-Optim_p_min)
Optim_index = (int)rint(Optim_bins * (p -Optim_p_min)/(Optim_p_max-Optim_p_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
if (Optim_Passing_p[Optim_index] >= Optim_Source_p[Optim_index])
{
Optim_Flag_Recycle |= OPTIM_MOD_P;
Optim_p += (Optim_index_p-Optim_index)*(Optim_p_max - Optim_p_min)/Optim_bins;
p = Optim_p;
}
else
Optim_index_p = Optim_index;
/* neutron is passed ! */
if (Optim_Source_vx[Optim_index_vx]
&& Optim_Source_vy[Optim_index_vy]
&& Optim_Source_vz[Optim_index_vz]
&& Optim_Source_x[Optim_index_x]
&& Optim_Source_y[Optim_index_y]
&& Optim_Source_s1[Optim_index_s1]
&& Optim_Source_s2[Optim_index_s2]
&& Optim_Source_p[Optim_index_p]
&& Optim_Reference_vx[Optim_index_vx]
&& Optim_Reference_vy[Optim_index_vy]
&& Optim_Reference_vz[Optim_index_vz]
&& Optim_Reference_x[Optim_index_x]
&& Optim_Reference_y[Optim_index_y]
&& Optim_Reference_s1[Optim_index_s1]
&& Optim_Reference_s2[Optim_index_s2]
&& Optim_Reference_p[Optim_index_p])
{
Optim_t1 = 1;
/* good neutrons have an improved distribution, so Ref/Source < 1 */
/* unmodified (form Ref kept fraction) neutrons have Passing < Ref*Optim_keep. their weight should be kept */
/* at the end there will be of those : 2*Optim_step*Optim_keep + (1-2*Optim_step)*Optim_keep = Optim_keep % of unmodified neutrons */
/* the remining part (1-Optim_keep neutrons) should have an integrated flux of (1-Optim_keep) */
Optim_t2 = (double)Optim_Reference_vx[Optim_index_vx]/Optim_Source_vx[Optim_index_vx];
if (Optim_t2 < 1) Optim_good_vx = Optim_index_vx;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_vy[Optim_index_vy]/Optim_Source_vy[Optim_index_vy];
if (Optim_t2 < 1) Optim_good_vy = Optim_index_vy;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_vz[Optim_index_vz]/Optim_Source_vz[Optim_index_vz];
if (Optim_t2 < 1) Optim_good_vz = Optim_index_vz;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_x[Optim_index_x]/Optim_Source_x[Optim_index_x];
if (Optim_t2 < 1) Optim_good_x = Optim_index_x;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_y[Optim_index_y]/Optim_Source_y[Optim_index_y];
if (Optim_t2 < 1) Optim_good_y = Optim_index_y;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_s1[Optim_index_s1]/Optim_Source_s1[Optim_index_s1];
if (Optim_t2 < 1) Optim_good_s1= Optim_index_s1;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_s2[Optim_index_s2]/Optim_Source_s2[Optim_index_s2];
if (Optim_t2 < 1) Optim_good_s2= Optim_index_s2;
Optim_t1 *= Optim_t2;
Optim_t2 = (double)Optim_Reference_p[Optim_index_p]/Optim_Source_p[Optim_index_p];
if (Optim_t2 < 1) Optim_good_p = Optim_index_p;
Optim_t1 *= Optim_t2;
/* now normalize to intial distribution */
Optim_p *= Optim_t1;
if (Optim_Flag_Recycle) { Optim_n_redirect++; } /* Optim_p /= Optim_Flag_Recycle; */
p = Optim_p;
}
else
ABSORB; /* can't modify neutron weight -> eject */
Optim_Passing_vx[Optim_index_vx]++;
Optim_Passing_vy[Optim_index_vy]++;
Optim_Passing_vz[Optim_index_vz]++;
Optim_Passing_x[Optim_index_x]++;
Optim_Passing_y[Optim_index_y]++;
Optim_Passing_s1[Optim_index_s1]++;
Optim_Passing_s2[Optim_index_s2]++;
Optim_Passing_p[Optim_index_p]++;
Optim_Passing_Counts++;
Optim_Passing_Flux += p;
} /* end if OPTIM_PHASE_OPTIM */
} /* end if xy in optimizer */
/* end trace */
%}
FINALLY
%{
if (strlen(file) > 0)
{
hfile = fopen(file, "w");
if(!hfile)
{
fprintf(stderr, "Error: %s : could not open output file '%s'\n", mccompcurname, file);
}
else
{
fprintf(hfile,"# Instrument-source: %s\n", mcinstrument_source);
mcruninfo_out("# ", hfile);
fprintf(hfile,"# type: array_2d(%i,6) \n",Optim_bins);
fprintf(hfile,"# component: %s\n", mccompcurname);
fprintf(hfile,"# title: General Optimizer distributions\n");
fprintf(hfile,"# filename: '%s'\n",file);
fprintf(hfile,"# variables: x dx y dy vx dvx vy dvy vz dvz s1 ds1 s2 ds2 p dp\n");
fprintf(hfile,"# xvar: (x y vx vy vz s1 s2 p)\n");
fprintf(hfile,"# yvar: (dx dy dvx dvy dvz ds1 ds2 dp)\n");
fprintf(hfile,"# xlabel: 'Distributions'\n");
fprintf(hfile,"# ylabel: 'Counts'\n");
if (Optim_Normal_Monitor_Counts != 0)
fprintf(hfile,"# Optimizer speedup estimate: %.3g [Monitor Normal counts %i (extrapolated), Optimized %i ]\n", (double)(Optim_Total_Monitor_Counts)/Optim_Normal_Monitor_Counts*2*Optim_step,(int)ceil(Optim_Normal_Monitor_Counts/2/Optim_step), Optim_Total_Monitor_Counts);
fprintf(hfile,"# Optimizer options: %s\n",options);
fprintf(hfile,"# Redirected neutrons: %i (%.2f \%)\n",Optim_n_redirect,(double)(100*Optim_n_redirect/mcget_ncount()));
fprintf(hfile,"# data: Optimzed Source\n");
for (Optim_index = 0; Optim_index < Optim_bins; Optim_index++)
{
fprintf(hfile,"%10.4g ",(Optim_x_min+((Optim_index+0.5)/Optim_bins)*(Optim_x_max - Optim_x_min)));
fprintf(hfile,"%10i\t",Optim_Source_x[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_y_min+((Optim_index+0.5)/Optim_bins)*(Optim_y_max - Optim_y_min)));
fprintf(hfile,"%10i\t",Optim_Source_y[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vx_min+((Optim_index+0.5)/Optim_bins)*(Optim_vx_max - Optim_vx_min)));
fprintf(hfile,"%10i\t",Optim_Source_vx[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vy_min+((Optim_index+0.5)/Optim_bins)*(Optim_vy_max - Optim_vy_min)));
fprintf(hfile,"%10i\t",Optim_Source_vy[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vz_min+((Optim_index+0.5)/Optim_bins)*(Optim_vz_max - Optim_vz_min)));
fprintf(hfile,"%10i\t",Optim_Source_vz[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s1_min+((Optim_index+0.5)/Optim_bins)*(Optim_s1_max - Optim_s1_min)));
fprintf(hfile,"%10i\t",Optim_Source_s1[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s2_min+((Optim_index+0.5)/Optim_bins)*(Optim_s2_max - Optim_s2_min)));
fprintf(hfile,"%10i\t",Optim_Source_s2[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_p_min+((Optim_index+0.5)/Optim_bins)*(Optim_p_max - Optim_p_min)));
fprintf(hfile,"%10i\t",Optim_Source_p[Optim_index]);
fprintf(hfile,"\n");
}
fprintf(hfile,"# data: Reference Source\n");
for (Optim_index = 0; Optim_index < Optim_bins; Optim_index++)
{
fprintf(hfile,"%10.4g ",(Optim_x_min+((Optim_index+0.5)/Optim_bins)*(Optim_x_max - Optim_x_min)));
fprintf(hfile,"%10i\t",Optim_Reference_x[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_y_min+((Optim_index+0.5)/Optim_bins)*(Optim_y_max - Optim_y_min)));
fprintf(hfile,"%10i\t",Optim_Reference_y[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vx_min+((Optim_index+0.5)/Optim_bins)*(Optim_vx_max - Optim_vx_min)));
fprintf(hfile,"%10i\t",Optim_Reference_vx[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vy_min+((Optim_index+0.5)/Optim_bins)*(Optim_vy_max - Optim_vy_min)));
fprintf(hfile,"%10i\t",Optim_Reference_vy[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vz_min+((Optim_index+0.5)/Optim_bins)*(Optim_vz_max - Optim_vz_min)));
fprintf(hfile,"%10i\t",Optim_Reference_vz[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s1_min+((Optim_index+0.5)/Optim_bins)*(Optim_s1_max - Optim_s1_min)));
fprintf(hfile,"%10i\t",Optim_Reference_s1[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s2_min+((Optim_index+0.5)/Optim_bins)*(Optim_s2_max - Optim_s2_min)));
fprintf(hfile,"%10i\t",Optim_Reference_s2[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_p_min+((Optim_index+0.5)/Optim_bins)*(Optim_p_max - Optim_p_min)));
fprintf(hfile,"%10i\t",Optim_Reference_p[Optim_index]);
fprintf(hfile,"\n");
}
fprintf(hfile,"# data: Passing\n");
for (Optim_index = 0; Optim_index < Optim_bins; Optim_index++)
{
fprintf(hfile,"%10.4g ",(Optim_x_min+((Optim_index+0.5)/Optim_bins)*(Optim_x_max - Optim_x_min)));
fprintf(hfile,"%10i\t",Optim_Passing_x[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_y_min+((Optim_index+0.5)/Optim_bins)*(Optim_y_max - Optim_y_min)));
fprintf(hfile,"%10i\t",Optim_Passing_y[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vx_min+((Optim_index+0.5)/Optim_bins)*(Optim_vx_max - Optim_vx_min)));
fprintf(hfile,"%10i\t",Optim_Passing_vx[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vy_min+((Optim_index+0.5)/Optim_bins)*(Optim_vy_max - Optim_vy_min)));
fprintf(hfile,"%10i\t",Optim_Passing_vy[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vz_min+((Optim_index+0.5)/Optim_bins)*(Optim_vz_max - Optim_vz_min)));
fprintf(hfile,"%10i\t",Optim_Passing_vz[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s1_min+((Optim_index+0.5)/Optim_bins)*(Optim_s1_max - Optim_s1_min)));
fprintf(hfile,"%10i\t",Optim_Passing_s1[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s2_min+((Optim_index+0.5)/Optim_bins)*(Optim_s2_max - Optim_s2_min)));
fprintf(hfile,"%10i\t",Optim_Passing_s2[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_p_min+((Optim_index+0.5)/Optim_bins)*(Optim_p_max - Optim_p_min)));
fprintf(hfile,"%10i\t",Optim_Passing_p[Optim_index]);
fprintf(hfile,"\n");
}
fprintf(hfile,"# data: New_Source\n");
for (Optim_index = 0; Optim_index < Optim_bins; Optim_index++)
{
fprintf(hfile,"%10.4g ",(Optim_x_min+((Optim_index+0.5)/Optim_bins)*(Optim_x_max - Optim_x_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_x[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_y_min+((Optim_index+0.5)/Optim_bins)*(Optim_y_max - Optim_y_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_y[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vx_min+((Optim_index+0.5)/Optim_bins)*(Optim_vx_max - Optim_vx_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_vx[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vy_min+((Optim_index+0.5)/Optim_bins)*(Optim_vy_max - Optim_vy_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_vy[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_vz_min+((Optim_index+0.5)/Optim_bins)*(Optim_vz_max - Optim_vz_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_vz[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s1_min+((Optim_index+0.5)/Optim_bins)*(Optim_s1_max - Optim_s1_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_s1[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_s2_min+((Optim_index+0.5)/Optim_bins)*(Optim_s2_max - Optim_s2_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_s2[Optim_index]);
fprintf(hfile,"%10.4g ",(Optim_p_min+((Optim_index+0.5)/Optim_bins)*(Optim_p_max - Optim_p_min)));
fprintf(hfile,"%10i\t",Optim_New_Source_p[Optim_index]);
fprintf(hfile,"\n");
}
fclose(hfile);
}
if (Optim_Flag_Verbose)
{
printf("End of optimization\n");
printf("Optim_Normal_Monitor_Counts = %i (2 steps), Optim_Total_Monitor_Counts = %i \n",Optim_Normal_Monitor_Counts, Optim_Total_Monitor_Counts);
if (Optim_Normal_Monitor_Counts != 0)
printf("Optimizer speedup : %.3g \n", (double)(Optim_Total_Monitor_Counts)/Optim_Normal_Monitor_Counts*2*Optim_step);
printf("Number of redirections : %i\n",Optim_n_redirect);
printf("Counts : reference = %i, passing = %i, monitor = %i\n", Optim_Reference_Counts, Optim_Passing_Counts, Optim_Monitor_Counts);
printf("Flux : reference = %.2g, passing = %.2g, monitor = %.2g\n", Optim_Reference_Flux, Optim_Passing_Flux, Optim_Monitor_Flux);
}
}
%}
MCDISPLAY
%{
magnify("xy");
multiline(5, (double)xmin, (double)ymin, 0.0,
(double)xmax, (double)ymin, 0.0,
(double)xmax, (double)ymax, 0.0,
(double)xmin, (double)ymax, 0.0,
(double)xmin, (double)ymin, 0.0);
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, version 1.1, released
* Maintained by Kristian Nielsen and Kim Lefmann,
* Risoe National Laboratory, Roskilde, Denmark
*
* Component: Monitor_Optimizer
*
* Written by: EF, 17 Sept 1999
*
* A component that optimizes the neutron flux passing through the
* Source_optimizer in order to have the maximum flux at the
* Monitor_Optimizer position.
* Source_optimizer should be placed just after the source.
* Monitor_Optimizer should be placed at position to optimize.
*
* INPUT PARAMETERS:
*
* xmin: Lower x bound of detector opening (m)
* xmax: Upper x bound of detector opening (m)
* ymin: Lower y bound of detector opening (m)
* ymax: Upper y bound of detector opening (m)
*
* OUTPUT PARAMETERS:
*
* distributions, optimizer efficiency
*
*******************************************************************************/
DEFINE COMPONENT Monitor_Optimizer
DEFINITION PARAMETERS (xmin, xmax, ymin, ymax)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
#ifndef MCSTAS_GEN_OPTIM
#error McStas : Source_Optimizer component has to be used before Monitor_Optimizer
#endif
%}
INITIALIZE
%{
%}
TRACE
%{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Optim_Monitor_Counts++; /* initialized in OPTIM_PHASE_SET_REF */
Optim_Monitor_Flux += p;
Optim_Total_Monitor_Counts++;
if (Optim_Flag_Auto
&& (Optim_Phase == OPTIM_PHASE_GET_REF))
{
if (((Optim_Reference_Counts > 2*Optim_Phase_Counts*Optim_step) && (Optim_Monitor_Counts >= Optim_bins*5))
|| (Optim_Monitor_Counts >= Optim_bins*10))
{
Optim_Phase_Counts_R = 0; /* enough counts on monitor */
if (Optim_Flag_Verbose)
{
printf(">> AUTO monitor has reached %i counts (non optimized",Optim_Monitor_Counts);
if (Optim_Flag_Absorb) printf(", absorbed");
printf(")\n");
}
Optim_step = (double)Optim_Reference_Counts/Optim_Phase_Counts;
Optim_Phase_Counts = Optim_Reference_Counts;
}
}
if ((Optim_Phase == OPTIM_PHASE_GET_REF)
|| ((Optim_Phase == OPTIM_PHASE_OPTIM) && (Optim_Flag_Continuous) )) /* build the Optimized Source distributions */
{
if (Optim_Monitor_pmax == 0 & p > Optim_Monitor_pmax) Optim_Monitor_pmax = p;
if (Optim_vx_max-Optim_vx_min)
Optim_index = (int)rint(Optim_bins * (Optim_vx -Optim_vx_min)/(Optim_vx_max-Optim_vx_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_vx[Optim_index]++;
if (Optim_vy_max-Optim_vy_min)
Optim_index = (int)rint(Optim_bins * (Optim_vy -Optim_vy_min)/(Optim_vy_max-Optim_vy_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_vy[Optim_index]++;
if (Optim_vz_max-Optim_vz_min)
Optim_index = (int)rint(Optim_bins * (Optim_vz -Optim_vz_min)/(Optim_vz_max-Optim_vz_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_vz[Optim_index]++;
if (Optim_x_max-Optim_x_min)
Optim_index = (int)rint(Optim_bins * (Optim_x -Optim_x_min)/(Optim_x_max-Optim_x_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_x[Optim_index]++;
if (Optim_y_max-Optim_y_min)
Optim_index = (int)rint(Optim_bins * (Optim_y -Optim_y_min)/(Optim_y_max-Optim_y_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_y[Optim_index]++;
if (Optim_s1_max-Optim_s1_min)
Optim_index = (int)rint(Optim_bins * (Optim_s1 -Optim_s1_min)/(Optim_s1_max-Optim_s1_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_s1[Optim_index]++;
if (Optim_s2_max-Optim_s2_min)
Optim_index = (int)rint(Optim_bins * (Optim_s2 -Optim_s2_min)/(Optim_s2_max-Optim_s2_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_s2[Optim_index]++;
if (Optim_p_max-Optim_p_min)
Optim_index = (int)rint(Optim_bins * (Optim_p -Optim_p_min)/(Optim_p_max-Optim_p_min));
else
Optim_index = 0;
if (Optim_index < 0) Optim_index = 0;
if (Optim_index >= Optim_bins) Optim_index = Optim_bins - 1;
Optim_New_Source_p[Optim_index]++;
} /* end if Optim_Phase */
if (Optim_Flag_Absorb && (Optim_Phase == OPTIM_PHASE_GET_REF))
ABSORB;
} /* end if xy in optimizer */
/* end trace */
%}
FINALLY
%{
/* initial Reference distribution arrays (for weights) */
free(Optim_Reference_x);
free(Optim_Reference_y);
free(Optim_Reference_vx);
free(Optim_Reference_vy);
free(Optim_Reference_vz);
free(Optim_Reference_s1);
free(Optim_Reference_s2);
free(Optim_Reference_p);
/* optimized Source distribution arrays (to reach) */
free(Optim_Source_x);
free(Optim_Source_y);
free(Optim_Source_vx);
free(Optim_Source_vy);
free(Optim_Source_vz);
free(Optim_Source_s1);
free(Optim_Source_s2);
free(Optim_Source_p);
/* optimized New_Source distribution arrays (to reach in next step, passed to Source) */
free(Optim_New_Source_x);
free(Optim_New_Source_y);
free(Optim_New_Source_vx);
free(Optim_New_Source_vy);
free(Optim_New_Source_vz);
free(Optim_New_Source_s1);
free(Optim_New_Source_s2);
free(Optim_New_Source_p);
/* Passing distribution arrays (should grow to reach Source) */
free(Optim_Passing_x);
free(Optim_Passing_y);
free(Optim_Passing_vx);
free(Optim_Passing_vy);
free(Optim_Passing_vz);
free(Optim_Passing_s1);
free(Optim_Passing_s2);
free(Optim_Passing_p);
%}
MCDISPLAY
%{
magnify("xy");
multiline(5, (double)xmin, (double)ymin, 0.0,
(double)xmax, (double)ymin, 0.0,
(double)xmax, (double)ymax, 0.0,
(double)xmin, (double)ymax, 0.0,
(double)xmin, (double)ymin, 0.0);
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, version 1.0, released October 26, 1998
* Maintained by Kristian Nielsen and Kim Lefmann,
* Risoe National Laboratory, Roskilde, Denmark
*
* Component: Guide2.
*
* Written by: KN, September 2 1998
* Modified by: KL, October 6, 1998
*
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* For details on the geometry calculation see the description in the McStas
* reference manual.
*
* INPUT PARAMETERS:
*
* w1: (m) Width at the guide entry
* h1: (m) Height at the guide entry
* w2: (m) Width at the guide exit
* h2: (m) Height at the guide exit
* l: (m) length of guide
* R0: (1) Low-angle reflectivity
* Qc: (AA-1) Critical scattering vector
* alpha: (AA) Slope of reflectivity
* mh: (1) m-horizontal value of material
* mv: (1) m-vertical value of material
* W: (AA-1) Width of supermirror cut-off
*
* Example values: m=4 Qc=0.02 W=1/300 alpha=6.49 R0=1
*/
DEFINE COMPONENT Guide2
DEFINITION PARAMETERS (w1, h1, w2, h2, l, R0, Qc, alpha, mh, mv, W)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
TRACE
%{
double t1,t2; /* Intersection times. */
double av,ah,bv,bh,cv1,cv2,ch1,ch2,d; /* Intermediate values */
double vdotn_v1,vdotn_v2,vdotn_h1,vdotn_h2; /* Dot products. */
int i; /* Which mirror hit? */
double q; /* Q [1/AA] of reflection */
double vlen2,nlen2; /* Vector lengths squared */
/* ToDo: These could be precalculated. */
double ww = .5*(w2 - w1), hh = .5*(h2 - h1);
double whalf = .5*w1, hhalf = .5*h1;
double lwhalf = l*whalf, lhhalf = l*hhalf;
/* Propagate neutron to guide entrance. */
PROP_Z0;
if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
ABSORB;
for(;;)
{
/* Compute the dot products of v and n for the four mirrors. */
av = l*vx; bv = ww*vz;
ah = l*vy; bh = hh*vz;
vdotn_v1 = bv + av; /* Left vertical */
vdotn_v2 = bv - av; /* Right vertical */
vdotn_h1 = bh + ah; /* Lower horizontal */
vdotn_h2 = bh - ah; /* Upper horizontal */
/* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */
cv1 = -whalf*l - z*ww; cv2 = x*l;
ch1 = -hhalf*l - z*hh; ch2 = y*l;
/* Compute intersection times. */
t1 = (l - z)/vz;
i = 0;
if(vdotn_v1 < 0 && (t2 = (cv1 - cv2)/vdotn_v1) < t1)
{
t1 = t2;
i = 1;
}
if(vdotn_v2 < 0 && (t2 = (cv1 + cv2)/vdotn_v2) < t1)
{
t1 = t2;
i = 2;
}
if(vdotn_h1 < 0 && (t2 = (ch1 - ch2)/vdotn_h1) < t1)
{
t1 = t2;
i = 3;
}
if(vdotn_h2 < 0 && (t2 = (ch1 + ch2)/vdotn_h2) < t1)
{
t1 = t2;
i = 4;
}
if(i == 0)
break; /* Neutron left guide. */
PROP_DT(t1);
switch(i)
{
case 1: /* Left vertical mirror */
nlen2 = l*l + ww*ww;
q = V2Q*(-2)*vdotn_v1/sqrt(nlen2);
d = 2*vdotn_v1/nlen2;
vx = vx - d*l;
vz = vz - d*ww;
break;
case 2: /* Right vertical mirror */
nlen2 = l*l + ww*ww;
q = V2Q*(-2)*vdotn_v2/sqrt(nlen2);
d = 2*vdotn_v2/nlen2;
vx = vx + d*l;
vz = vz - d*ww;
break;
case 3: /* Lower horizontal mirror */
nlen2 = l*l + hh*hh;
q = V2Q*(-2)*vdotn_h1/sqrt(nlen2);
d = 2*vdotn_h1/nlen2;
vy = vy - d*l;
vz = vz - d*hh;
break;
case 4: /* Upper horizontal mirror */
nlen2 = l*l + hh*hh;
q = V2Q*(-2)*vdotn_h2/sqrt(nlen2);
d = 2*vdotn_h2/nlen2;
vy = vy + d*l;
vz = vz - d*hh;
break;
}
/* Now compute reflectivity. */
if(q > Qc)
{
double arg = (q - (i <= 2 ? mv : mh)*Qc)/W;
if(arg < 10)
p *= .5*(1-tanh(arg))*(1-alpha*(q-Qc));
else
ABSORB; /* Cutoff ~ 1E-10 */
}
p *= R0;
}
%}
END
-------------- next part --------------
DEFINE INSTRUMENT IN14(KI,WN,ORDER,MHV)
/* preprocess with : mcstas in14_4.instr */
/* compile on mica with : cc -Ofast -64 -o in14_3 in14_3.c -lm */
/* Test with : */
/* in14_6 -n 1e6 KI=2.66 WN=0.03 ORDER=1 MHV=3 */
DECLARE
%{
#define VERSION "60"
double L1 = 16.68; /* source-mono */
double L2 = 2.12; /* mono-sample */
double L3 = 1.35; /* sample-ana */
double L4 = 0.70; /* ana-det */
int SM,SS,SA;
double A1,A3,A5;
double LM1, LM1b; /* distances to monitors M1 and M2 */
double LM2, LM2b;
double A2,A4,A6,RM,RH;
char *pfile;
char *efile;
char *dfile;
char *kfile;
/* ==================== Source description ==================== */
double EN;
double D_EN;
double EMIN, EMAX;
double FLUX = 1e13; /* n/cm^2/s on guide entry */
/* ==================== Monochromator desciption ==================== */
double ETAM = 30; /* Mosaicity used on monochromator 30 */
double DM = 3.355; /* Monochromator d-spacing in Angs */
/* PG002 Orders : 1st 3.355 2e 1.6775, 3e 1.1183 */
double mono_r0 = 0.9; /* mean reflectivity */
double mono_width = 0.15;
double mono_heigh = .122;
double mono_gap = 0; /* slates are adjacent */
int mono_segment_number_vert = 7;
int mono_segment_number_horz = 1;
double mono_curv_vert; /* Vertical Rotation between adjacent slabs. */
double mono_curv_horz; /* Horizontal Rotation between adjacent slabs. */
double mono_slate_heigh; /* size (height) of slates */
double mono_slate_width; /* size (width) of slates */
double mono_q; /* Q vector for bragg scattering with monochromator and analysator */
double Ki, Ei;
double TM, GM;
/* ==================== Sample desciption ==================== */
double TU, TL;
double GU, GL;
/* ==================== Analyser desciption ==================== */
double ETAA = 30; /* Mosaicity used on analyser 30 */
double DA = 3.355; /* analyser d-spacing in Angs */
/* PG002 Orders : 1st 3.355 2e 1.6775, 3e 1.1183 */
double ana_r0 = 0.9; /* mean reflectivity */
double ana_width = 0.10;
double ana_heigh = .14;
double ana_gap = 0; /* slates are adjacent */
int ana_segment_number_vert = 7;
int ana_segment_number_horz = 1;
double ana_curv_vert; /* Vertical Rotation between adjacent slabs. */
double ana_curv_horz; /* Horizontal Rotation between adjacent slabs. */
double ana_slate_heigh; /* size (height) of slates */
double ana_slate_width; /* size (width) of slates */
double ana_q; /* Q vector for bragg scattering with monochromator and analysator */
double Kf, Ef;
double TA, GA;
%}
INITIALIZE
%{
double vi,vf,sample_q;
mono_q = 2*PI*ORDER/DM; /* Q mono in Angs-1 */
A4 = 0;
A2 = asin(mono_q/2/KI)*RAD2DEG*2;
A6 = A2;
printf("Instrument : IN14, v%s (21/09/99) on %s.\n",VERSION,getenv("HOSTNAME"));
/* SM : scattering at mono to the right (-1)/left(+1) */
/* SS : scattering at sample to the right (-1)/left(+1) */
/* SA : scattering at analyser to the right (-1)/left(+1) */
SM = 1; SS = 1; SA = 1;
A2 *= SM; /* A1 : mono theta (crystal) */
A1 = A2/2; /* A2 : mono 2 theta (arm to sample) */
A4 *= SS; /* A3 : sample theta */
A3 = A4/2; /* A4 : sample 2 theta (arm to analyser) */
A6 *= SA; /* A5 : analyser theta (crystal) */
A5 = A6/2; /* A6 : analyser 2 theta (arm to detector) */
TM = 0; /* TM : translation mono */
GM = 0; /* GM : tilt mono */
GU = 0; /* GU : tilt sample Up */
GL = 0; /* GL : tilt sample Low */
TU = 0; /* TU : translation sample Up */
TL = 0; /* TL : translation sample Low */
TA = 0; /* TA : translation analyser */
GA = 0; /* GA : tilt analyser */
/* RA : horizontal curvature analyser */
if ((fabs(mono_q/2/KI) < 1) && (sin(DEG2RAD*A1) != 0))
Ki = mono_q/2/sin(DEG2RAD*A1);
else
{
printf("Warn : Can't define incident wave vector Ki\n");
Ki = 0;
printf("Skipping simulation\n");
exit(-1);
}
vi = K2V*fabs(Ki);
Ei = VS2E*vi*vi;
EN = Ei; D_EN = .5; /* optimize source on Ki */
ana_q = 2*PI/DA; /* Q ana in Angs-1 */
if (sin(DEG2RAD*A5) != 0)
Kf = ana_q/2/sin(DEG2RAD*A5);
else
{
printf("Warn : Can't define outgoing wave vector Kf\n");
Kf = 0;
}
vf = K2V*fabs(Kf);
Ef = VS2E*vf*vf;
sample_q = sqrt(Ki*Ki + Kf*Kf -2*fabs(Ki*Kf)*cos(DEG2RAD*A4));
mono_slate_heigh = mono_heigh/mono_segment_number_vert; /* slates are adjacent */
mono_curv_vert = fabs(mono_slate_heigh*RAD2DEG/(2*L2*sin(DEG2RAD*A1))); /* RM : vertical mono curvature */
mono_slate_width = mono_width/mono_segment_number_horz; /* slates are adjacent */
mono_curv_horz = fabs(mono_slate_width*RAD2DEG*sin(DEG2RAD*A1)/(2*L2)); /* RH : horizontal mono curvature */
ana_slate_heigh = ana_heigh/ana_segment_number_vert; /* slates are adjacent */
ana_curv_vert = fabs(ana_slate_heigh*RAD2DEG/(2*L3*sin(DEG2RAD*A5))); /* RA : vertical ana curvature */
ana_slate_width = ana_width/ana_segment_number_horz; /* slates are adjacent */
ana_curv_horz = fabs(ana_slate_width*RAD2DEG*sin(DEG2RAD*A5)/(2*L3)); /* RHA : horizontal ana curvature */
/* print instrument config */
printf("Flat source, m=%.2f noze, width %.2f\n",MHV,WN);
printf("Monochromator : (DM = %g)\n",DM);
printf("A1 = %.2f, A2 = %.2f (deg)\n",A1,A2);
printf("Ki = %.3g Angs-1 (Energy = %.3g meV, Velocity = %.3g m/s) \n", Ki, Ei,vi);
printf("RM = %.3g Deg, RH = %.3g Deg\n",mono_curv_vert,mono_curv_horz);
printf("\n");
printf("Sample :\n");
printf("A3 = %.2f, A4 = %.2f (deg)\n",A3,A4);
printf("Energy transfert %.3g meV, Moment transfert %.3g Angs-1\n",Ef-Ei,sample_q);
printf("\n");
printf("Analyser : (DA = %g)\n",DA);
printf("A5 = %.2f, A6 = %.2f (deg)\n",A5,A6);
printf("Kf = %.3g Angs-1 (Energy = %.3g meV, Velocity = %.3g m/s) \n", Kf, Ef,vf);
printf("RA = %.3g Deg\n",ana_curv_vert);
printf("\n");
printf("Detectors :\n");
/* local variables ------------------------------------ */
LM1 = L2*.9; LM1b = LM1+0.001;
LM2 = L3/2; LM2b = LM2+0.001;
EMIN = EN - D_EN;
EMAX = EN + D_EN;
pfile = (char *)malloc(256);
efile = (char *)malloc(256);
dfile = (char *)malloc(256);
kfile = (char *)malloc(256);
sprintf(pfile,"sim/i%s_k%iw%id%im%i.psd",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
sprintf(efile,"sim/i%s_k%iw%id%im%i.nrj",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
sprintf(dfile,"sim/i%s_k%iw%id%im%i.div",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
sprintf(kfile,"sim/i%s_k%iw%id%im%i.kw",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
%}
TRACE
COMPONENT origin = Arm()
AT (0,0,0) ABSOLUTE
COMPONENT source = Source_flux(
radius = 0.20,
dist = 2.16,
xw = 0.06, yh = 0.12,
E0 = EN,
dE = D_EN,
flux=FLUX)
AT (0,0,0) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
COMPONENT optim_s = Source_Optimizer(
xmin = -0.03, xmax = 0.03,
ymin = -0.06, ymax = 0.06,
bins = -1, step = 10, keep = 10,
file="source.optim",
options="verbose+absorb+auto")
AT (0,0,2.14) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
COMPONENT kw1 = kw_monitor(
xmin = -0.01, xmax = 0.01,
ymin = -0.01, ymax = 0.01,
filename = "kguide.kw")
AT(0, 0, 2.145) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
COMPONENT doigt_de_gant = Guide(
w1 = 0.06, h1 = 0.12,
w2 = 0.06, h2 = 0.12,
l = 2.75, /* guide into the doigt de gant */
R0 = 1, m=1.2, /* Ni 58 */
Qc = 0.021, alpha = 2.33, W = 2e-3)
AT (0,0,2.15) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
COMPONENT external_guide = Guide2(
w1 = 0.06, h1 = 0.12,
w2 = 0.06, h2 = 0.12,
l = 13.67, /* external guide between doigt de gant and mono */
R0 = 1,
Qc = 0.021, alpha = 2.33,
mh=1.2, /* 1.2 Ni 58, 3 Super mirror */
mv=1.2,
W = 2e-3)
AT (0,0,4.91) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
/* -------------- Start of monochromator building -------------- */
/* support of mono */
COMPONENT focus_mono = Arm()
AT (0, 0, L1) RELATIVE origin ROTATED (0, A1, 0) RELATIVE origin
COMPONENT m0 = Mon_2foc(
zwidth=mono_slate_width,
ywidth=mono_slate_heigh,
gap=mono_gap,
NH=mono_segment_number_vert,
NV=mono_segment_number_horz,
mosaich=ETAM,
mosaicv=ETAM,
r0=mono_r0, Q=mono_q,
RH=mono_curv_vert,
RV=mono_curv_horz)
AT (TM, 0, 0) RELATIVE focus_mono
ROTATED (0, 0, GM) RELATIVE focus_mono
/* on mono, pointing towards sample */
COMPONENT out_mono = Arm()
AT (0,0,0) RELATIVE focus_mono ROTATED (0, A2, 0) RELATIVE origin
/* -------------- End of monochromator building -------------- */
COMPONENT noze = Guide2(
w1 = 0.05, h1 = 0.05,
w2 = WN, h2 = 0.05,
l = .825,
R0 = 1,
Qc = 0.021, alpha = 2.33,
mh=1e-5, mv = MHV, /* Ni 58 */
W = 2e-3)
AT (0, 0, .9) RELATIVE out_mono ROTATED (0,0,0) RELATIVE out_mono
/*
COMPONENT M1 = Monitor(
xmin = -0.03, xmax = 0.03,
ymin = -0.06, ymax = 0.06)
AT (0, 0, LM1) RELATIVE out_mono ROTATED (0,0,0) RELATIVE out_mono
*/
/* -------------- Start of sample building -------------- */
/* support of sample */
COMPONENT focus_sample = Arm()
AT (0, 0, L2) RELATIVE out_mono ROTATED (0,A3,0) RELATIVE out_mono
/*
COMPONENT sample = Powder1(
radius = 0.007,
h = 0.015,
q = 1.8049,
d_phi0 = 4,
pack = 1, j = 6, DW = 1,
F2 = 56.8,
Vc = 85.0054, sigma_a = 0.463,
target_x = alu_focus_x,
target_y = 0, target_z = 1000)
AT (GU, 0, GL) RELATIVE focus_sample ROTATED (TU,0,TL) RELATIVE focus_sample
*/
/*
COMPONENT AtSample = Monitor(
xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05)
AT (0, 0, 0) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
*/
COMPONENT optim_m = Monitor_Optimizer(
xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05)
AT(0, 0, 0) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT PSDSample = PSD_monitor(
xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05,
nx = 50, ny = 50,
filename = pfile)
AT(0, 0, 0.001) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT Int4Sample = Monitor(
xmin = -0.02, xmax = 0.02,
ymin = -0.02, ymax = 0.02)
AT(0, 0, 0.002) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT Flux1Sample = Monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005)
AT(0, 0, 0.003) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT ESample = E_monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005,
Emin = EMIN, Emax = EMAX, nchan = 21,
filename = efile)
AT(0, 0, 0.004) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT DivSample = Divergence_monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005,
nv = 10, nh= 10,
v_maxdiv = 1, h_maxdiv = 1,
filename = dfile)
AT(0, 0, 0.005) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT kw = kw_monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005,
filename = kfile)
AT(0, 0, 0.006) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
/* on sample, pointing towards ana */
COMPONENT out_sample = Arm()
AT (0,0,0) RELATIVE focus_sample ROTATED (0, A4, 0) RELATIVE out_mono
/* -------------- End of sample building -------------- */
/*
COMPONENT M2 = Monitor(
xmin = -0.1, xmax = 0.1,
ymin = -0.1, ymax = 0.1)
AT (0, 0, LM2) RELATIVE out_sample ROTATED (0,0,0) RELATIVE out_sample
*/
/* -------------- Start of analyzer building -------------- */
/* support of analyzer */
/*
COMPONENT focus_ana = Arm()
AT (0, 0, L3) RELATIVE out_sample ROTATED (0,A5,0) RELATIVE out_sample
COMPONENT a0 = Mon_2foc(
zwidth = ana_half_width,
ywidth = ana_half_heigh,
gap = 0,
NH = 1,
NV = 1,
mosaich = ETAA,
mosaicv = ETAA,
r0 = ana_r0,
Q = ana_q,
RH=ana_curv_vert,
RV=ana_curv_horz)
AT (TA, 0, 0) RELATIVE focus_ana ROTATED (0, 0, GA) RELATIVE focus_ana
COMPONENT out_ana = Arm()
AT (0, 0, 0) RELATIVE focus_ana ROTATED (0, A6, 0) RELATIVE out_sample
*/
/* -------------- End of analyzer building -------------- */
/*
COMPONENT focus_det = Arm()
AT (0, 0, L4) RELATIVE out_ana ROTATED (0,0,0) RELATIVE out_ana
COMPONENT Detector = Monitor(
xmin = -0.02, xmax = 0.02,
ymin = -0.05, ymax = 0.05)
AT(0, 0, 0) RELATIVE focus_det ROTATED (0,0,0) RELATIVE focus_det
COMPONENT emon2 = E_monitor(
xmin = -0.02, xmax = 0.02,
ymin = -0.05, ymax = 0.05,
Emin = 10, Emax = 20, nchan = 35,
filename = "sim/in14_EM2.vmon")
AT(0, 0, 0.01) RELATIVE focus_det ROTATED (0,0,0) RELATIVE focus_det
*/
FINALLY
%{
printf("Output files : %s %s %s %s\n",pfile,efile,dfile,kfile);
free(pfile);
free(dfile);
free(efile);
free(kfile);
%}
END
-------------- next part --------------
DEFINE INSTRUMENT IN14(KI,WN,ORDER,MHV)
/* preprocess with : mcstas in14_4.instr */
/* compile on mica with : cc -Ofast -64 -o in14_3 in14_3.c -lm */
/* Test with : in14_4 KI=1.55 WN=5 */
/* in14_6n -n 1e6 KI=2.66 WN=0.03 ORDER=1 MHV=3 */
DECLARE
%{
#define VERSION "61"
double L1 = 16.68; /* source-mono */
double L2 = 2.12; /* mono-sample */
double L3 = 1.35; /* sample-ana */
double L4 = 0.70; /* ana-det */
int SM,SS,SA;
double A1,A3,A5;
double LM1, LM1b; /* distances to monitors M1 and M2 */
double LM2, LM2b;
double A2,A4,A6,RM,RH;
char *pfile;
char *efile;
char *dfile;
char *kfile;
/* ==================== Source description ==================== */
double EN;
double D_EN;
double EMIN, EMAX;
double FLUX = 1e13; /* n/cm^2/s on guide entry */
/* ==================== Monochromator desciption ==================== */
double ETAM = 30; /* Mosaicity used on monochromator 30 */
double DM = 3.355; /* Monochromator d-spacing in Angs */
/* PG002 Orders : 1st 3.355 2e 1.6775, 3e 1.1183 */
double mono_r0 = 0.9; /* mean reflectivity */
double mono_width = 0.15;
double mono_heigh = .122;
double mono_gap = 0; /* slates are adjacent */
int mono_segment_number_vert = 7;
int mono_segment_number_horz = 1;
double mono_curv_vert; /* Vertical Rotation between adjacent slabs. */
double mono_curv_horz; /* Horizontal Rotation between adjacent slabs. */
double mono_slate_heigh; /* size (height) of slates */
double mono_slate_width; /* size (width) of slates */
double mono_q; /* Q vector for bragg scattering with monochromator and analysator */
double Ki, Ei;
double TM, GM;
/* ==================== Sample desciption ==================== */
double TU, TL;
double GU, GL;
/* ==================== Analyser desciption ==================== */
double ETAA = 30; /* Mosaicity used on analyser 30 */
double DA = 3.355; /* analyser d-spacing in Angs */
/* PG002 Orders : 1st 3.355 2e 1.6775, 3e 1.1183 */
double ana_r0 = 0.9; /* mean reflectivity */
double ana_width = 0.10;
double ana_heigh = .14;
double ana_gap = 0; /* slates are adjacent */
int ana_segment_number_vert = 7;
int ana_segment_number_horz = 1;
double ana_curv_vert; /* Vertical Rotation between adjacent slabs. */
double ana_curv_horz; /* Horizontal Rotation between adjacent slabs. */
double ana_slate_heigh; /* size (height) of slates */
double ana_slate_width; /* size (width) of slates */
double ana_q; /* Q vector for bragg scattering with monochromator and analysator */
double Kf, Ef;
double TA, GA;
%}
INITIALIZE
%{
double vi,vf,sample_q;
mono_q = 2*PI*ORDER/DM; /* Q mono in Angs-1 */
A4 = 0;
A2 = asin(mono_q/2/KI)*RAD2DEG*2;
A6 = A2;
printf("Instrument : IN14, v%s (21/09/99) on %s.\n",VERSION,getenv("HOSTNAME"));
/* SM : scattering at mono to the right (-1)/left(+1) */
/* SS : scattering at sample to the right (-1)/left(+1) */
/* SA : scattering at analyser to the right (-1)/left(+1) */
SM = 1; SS = 1; SA = 1;
A2 *= SM; /* A1 : mono theta (crystal) */
A1 = A2/2; /* A2 : mono 2 theta (arm to sample) */
A4 *= SS; /* A3 : sample theta */
A3 = A4/2; /* A4 : sample 2 theta (arm to analyser) */
A6 *= SA; /* A5 : analyser theta (crystal) */
A5 = A6/2; /* A6 : analyser 2 theta (arm to detector) */
TM = 0; /* TM : translation mono */
GM = 0; /* GM : tilt mono */
GU = 0; /* GU : tilt sample Up */
GL = 0; /* GL : tilt sample Low */
TU = 0; /* TU : translation sample Up */
TL = 0; /* TL : translation sample Low */
TA = 0; /* TA : translation analyser */
GA = 0; /* GA : tilt analyser */
/* RA : horizontal curvature analyser */
if ((fabs(mono_q/2/KI) < 1) && (sin(DEG2RAD*A1) != 0))
Ki = mono_q/2/sin(DEG2RAD*A1);
else
{
printf("Warn : Can't define incident wave vector Ki\n");
Ki = 0;
printf("Skipping simulation\n");
exit(-1);
}
vi = K2V*fabs(Ki);
Ei = VS2E*vi*vi;
EN = Ei; D_EN = .5; /* optimize source on Ki */
ana_q = 2*PI/DA; /* Q ana in Angs-1 */
if (sin(DEG2RAD*A5) != 0)
Kf = ana_q/2/sin(DEG2RAD*A5);
else
{
printf("Warn : Can't define outgoing wave vector Kf\n");
Kf = 0;
}
vf = K2V*fabs(Kf);
Ef = VS2E*vf*vf;
sample_q = sqrt(Ki*Ki + Kf*Kf -2*fabs(Ki*Kf)*cos(DEG2RAD*A4));
mono_slate_heigh = mono_heigh/mono_segment_number_vert; /* slates are adjacent */
mono_curv_vert = fabs(mono_slate_heigh*RAD2DEG/(2*L2*sin(DEG2RAD*A1))); /* RM : vertical mono curvature */
mono_slate_width = mono_width/mono_segment_number_horz; /* slates are adjacent */
mono_curv_horz = fabs(mono_slate_width*RAD2DEG*sin(DEG2RAD*A1)/(2*L2)); /* RH : horizontal mono curvature */
ana_slate_heigh = ana_heigh/ana_segment_number_vert; /* slates are adjacent */
ana_curv_vert = fabs(ana_slate_heigh*RAD2DEG/(2*L3*sin(DEG2RAD*A5))); /* RA : vertical ana curvature */
ana_slate_width = ana_width/ana_segment_number_horz; /* slates are adjacent */
ana_curv_horz = fabs(ana_slate_width*RAD2DEG*sin(DEG2RAD*A5)/(2*L3)); /* RHA : horizontal ana curvature */
/* print instrument config */
printf("Flat source, m=%.2f noze, width %.2f\n",MHV,WN);
printf("Monochromator : (DM = %g)\n",DM);
printf("A1 = %.2f, A2 = %.2f (deg)\n",A1,A2);
printf("Ki = %.3g Angs-1 (Energy = %.3g meV, Velocity = %.3g m/s) \n", Ki, Ei,vi);
printf("RM = %.3g Deg, RH = %.3g Deg\n",mono_curv_vert,mono_curv_horz);
printf("\n");
printf("Sample :\n");
printf("A3 = %.2f, A4 = %.2f (deg)\n",A3,A4);
printf("Energy transfert %.3g meV, Moment transfert %.3g Angs-1\n",Ef-Ei,sample_q);
printf("\n");
printf("Analyser : (DA = %g)\n",DA);
printf("A5 = %.2f, A6 = %.2f (deg)\n",A5,A6);
printf("Kf = %.3g Angs-1 (Energy = %.3g meV, Velocity = %.3g m/s) \n", Kf, Ef,vf);
printf("RA = %.3g Deg\n",ana_curv_vert);
printf("\n");
printf("Detectors :\n");
/* local variables ------------------------------------ */
LM1 = L2*.9; LM1b = LM1+0.001;
LM2 = L3/2; LM2b = LM2+0.001;
EMIN = EN - D_EN;
EMAX = EN + D_EN;
pfile = (char *)malloc(256);
efile = (char *)malloc(256);
dfile = (char *)malloc(256);
kfile = (char *)malloc(256);
sprintf(pfile,"sim/i%s_k%iw%id%im%i.psd",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
sprintf(efile,"sim/i%s_k%iw%id%im%i.nrj",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
sprintf(dfile,"sim/i%s_k%iw%id%im%i.div",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
sprintf(kfile,"sim/i%s_k%iw%id%im%i.kw",VERSION,(int)floor(10*Ki+0.5),(int)floor(1000*WN),(int)floor(ORDER),(int)floor(MHV+0.5));
%}
TRACE
COMPONENT origin = Arm()
AT (0,0,0) ABSOLUTE
COMPONENT source = Source_flux(
radius = 0.20,
dist = 2.16,
xw = 0.06, yh = 0.12,
E0 = EN,
dE = D_EN,
flux=FLUX)
AT (0,0,0) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
COMPONENT doigt_de_gant = Guide(
w1 = 0.06, h1 = 0.12,
w2 = 0.06, h2 = 0.12,
l = 2.75, /* guide into the doigt de gant */
R0 = 1, m=1.2, /* Ni 58 */
Qc = 0.021, alpha = 2.33, W = 2e-3)
AT (0,0,2.15) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
COMPONENT external_guide = Guide2(
w1 = 0.06, h1 = 0.12,
w2 = 0.06, h2 = 0.12,
l = 13.67, /* external guide between doigt de gant and mono */
R0 = 1,
Qc = 0.021, alpha = 2.33,
mh=1.2, /* 1.2 Ni 58, 3 Super mirror */
mv=1.2,
W = 2e-3)
AT (0,0,4.91) RELATIVE origin ROTATED (0,0,0) RELATIVE origin
/* -------------- Start of monochromator building -------------- */
/* support of mono */
COMPONENT focus_mono = Arm()
AT (0, 0, L1) RELATIVE origin ROTATED (0, A1, 0) RELATIVE origin
COMPONENT m0 = Mon_2foc(
zwidth=mono_slate_width,
ywidth=mono_slate_heigh,
gap=mono_gap,
NH=mono_segment_number_vert,
NV=mono_segment_number_horz,
mosaich=ETAM,
mosaicv=ETAM,
r0=mono_r0, Q=mono_q,
RH=mono_curv_vert,
RV=mono_curv_horz)
AT (TM, 0, 0) RELATIVE focus_mono
ROTATED (0, 0, GM) RELATIVE focus_mono
/* on mono, pointing towards sample */
COMPONENT out_mono = Arm()
AT (0,0,0) RELATIVE focus_mono ROTATED (0, A2, 0) RELATIVE origin
/* -------------- End of monochromator building -------------- */
COMPONENT noze = Guide2(
w1 = 0.05, h1 = 0.05,
w2 = WN, h2 = 0.05,
l = .825,
R0 = 1,
Qc = 0.021, alpha = 2.33,
mh=1e-5, mv = MHV, /* Ni 58 */
W = 2e-3)
AT (0, 0, .9) RELATIVE out_mono ROTATED (0,0,0) RELATIVE out_mono
/*
COMPONENT M1 = Monitor(
xmin = -0.03, xmax = 0.03,
ymin = -0.06, ymax = 0.06)
AT (0, 0, LM1) RELATIVE out_mono ROTATED (0,0,0) RELATIVE out_mono
*/
/* -------------- Start of sample building -------------- */
/* support of sample */
COMPONENT focus_sample = Arm()
AT (0, 0, L2) RELATIVE out_mono ROTATED (0,A3,0) RELATIVE out_mono
/*
COMPONENT sample = Powder1(
radius = 0.007,
h = 0.015,
q = 1.8049,
d_phi0 = 4,
pack = 1, j = 6, DW = 1,
F2 = 56.8,
Vc = 85.0054, sigma_a = 0.463,
target_x = alu_focus_x,
target_y = 0, target_z = 1000)
AT (GU, 0, GL) RELATIVE focus_sample ROTATED (TU,0,TL) RELATIVE focus_sample
*/
/*
COMPONENT AtSample = Monitor(
xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05)
AT (0, 0, 0) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
*/
COMPONENT PSDSample = PSD_monitor(
xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05,
nx = 50, ny = 50,
filename = pfile)
AT(0, 0, 0.001) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT Int4Sample = Monitor(
xmin = -0.02, xmax = 0.02,
ymin = -0.02, ymax = 0.02)
AT(0, 0, 0.002) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT Flux1Sample = Monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005)
AT(0, 0, 0.003) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT ESample = E_monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005,
Emin = EMIN, Emax = EMAX, nchan = 21,
filename = efile)
AT(0, 0, 0.004) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT DivSample = Divergence_monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005,
nv = 10, nh= 10,
v_maxdiv = 1, h_maxdiv = 1,
filename = dfile)
AT(0, 0, 0.005) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
COMPONENT kw = kw_monitor(
xmin = -0.005, xmax = 0.005,
ymin = -0.005, ymax = 0.005,
filename = kfile)
AT(0, 0, 0.006) RELATIVE focus_sample ROTATED (0,0,0) RELATIVE focus_sample
/* on sample, pointing towards ana */
COMPONENT out_sample = Arm()
AT (0,0,0) RELATIVE focus_sample ROTATED (0, A4, 0) RELATIVE out_mono
/* -------------- End of sample building -------------- */
/*
COMPONENT M2 = Monitor(
xmin = -0.1, xmax = 0.1,
ymin = -0.1, ymax = 0.1)
AT (0, 0, LM2) RELATIVE out_sample ROTATED (0,0,0) RELATIVE out_sample
*/
/* -------------- Start of analyzer building -------------- */
/* support of analyzer */
/*
COMPONENT focus_ana = Arm()
AT (0, 0, L3) RELATIVE out_sample ROTATED (0,A5,0) RELATIVE out_sample
COMPONENT a0 = Mon_2foc(
zwidth = ana_half_width,
ywidth = ana_half_heigh,
gap = 0,
NH = 1,
NV = 1,
mosaich = ETAA,
mosaicv = ETAA,
r0 = ana_r0,
Q = ana_q,
RH=ana_curv_vert,
RV=ana_curv_horz)
AT (TA, 0, 0) RELATIVE focus_ana ROTATED (0, 0, GA) RELATIVE focus_ana
COMPONENT out_ana = Arm()
AT (0, 0, 0) RELATIVE focus_ana ROTATED (0, A6, 0) RELATIVE out_sample
*/
/* -------------- End of analyzer building -------------- */
/*
COMPONENT focus_det = Arm()
AT (0, 0, L4) RELATIVE out_ana ROTATED (0,0,0) RELATIVE out_ana
COMPONENT Detector = Monitor(
xmin = -0.02, xmax = 0.02,
ymin = -0.05, ymax = 0.05)
AT(0, 0, 0) RELATIVE focus_det ROTATED (0,0,0) RELATIVE focus_det
COMPONENT emon2 = E_monitor(
xmin = -0.02, xmax = 0.02,
ymin = -0.05, ymax = 0.05,
Emin = 10, Emax = 20, nchan = 35,
filename = "sim/in14_EM2.vmon")
AT(0, 0, 0.01) RELATIVE focus_det ROTATED (0,0,0) RELATIVE focus_det
*/
FINALLY
%{
printf("Output files : %s %s %s %s\n",pfile,efile,dfile,kfile);
free(pfile);
free(dfile);
free(efile);
%}
END
More information about the mcstas-users
mailing list