McStas : flux adapter
Farhi
farhi at ill.fr
Thu Oct 21 18:58:06 CEST 1999
Hello Kristian,
I've written a component in order to adapt a flux to a real reactor
distribution. It can handle tables versus energy, wavevector and
wavelength, and works with an external 'flux description' file in free
format (just needs to see two columns somewhere inside).
I attach an example (with the in14_6 instrument version "60a", also
provided). The optimizer works fine with the flux adapter . I've
improved the optimzer (works with any flux distrubution, auto mode is
now better, and some new options and default mode is more general). Hope
you'll enjoy this.
.In real usage, the flux description file should be normalized when in
'multiply' mode.
Concerning some 'spikes' of weight 'p' that I had after the optimization
in the (k,w) distribution, they are normal : as some neutrons are not
very efficient, and so not often used, thiri weight is bigger. One
should then integrate somehow the p distribution. A neutron alone is not
representative, but an assembly is ok.
Cheers... EF.
--
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/19991021/68f2403a/attachment.html>
-------------- 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 "60a"
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 fa = Flux_adapter(
xmin = -0.1, xmax = 0.1,
ymin = -0.1, ymax = 0.1,
file="source.flux",
options="")
AT (0,0,0.1) 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.14) 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="auto")
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 --------------
# energy multiply
0 3
2 3
10 3
100 3
-------------- next part --------------
/***********************************************************************
*
* McStas, version 1.1, released ??
* Maintained by Kristian Nielsen and Kim Lefmann,
* Risoe National Laboratory, Roskilde, Denmark
*
* Component: Flux_adapter
*
* Written by: EF, Oct 14, 1999
*
* The routine changes the neutron flux (weight) in order to match
* a reference source table on disk. This file can be on format
* (k[Angs-1],p), (omega[meV],p), (lambda[Angs],p) where p is the weight
* A linear interpolation is performed during simulation.
* This component should be placed after a source, in order to
* simulate a real source.
*
* file : name of the file to look at (two columns data)
* data should be sorted (ascending order) and monotonic
* file can contain options (see below) as comment
* options : string that can contain
* "[ k p ]" or "wavector" for file type
* "[ omega p]" or "energy"
* "[ lambda p ]" or "wavelength"
* "set" to set the weight according to the table
* "multiply" to multiply (instead of set) the weight by factor
* "add" to add to current flux
* "unactivate" to unactivate flux_adapter (in flux file for test)
* "verbose" to display additional informations
*
* EXAMPLE :
* (xmin = -0.1, xmax = 0.1,
* ymin = -0.1, ymax = 0.1,
* file="source.flux",
* options="")
* With file "source.flux" beeing something like
* # energy multiply
* # [ meV flux_factor ]
* 0 1
* 2 1.2
* 10 1.5
* 100 1
*
***********************************************************************/
DEFINE COMPONENT Flux_adapter
DEFINITION PARAMETERS (xmin, xmax, ymin, ymax, file, options)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (KE_Table,Weight_Table,Type_table,Mode_Table,Length_Table)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
#define LINE_MAX_LENGTH 1024
#define UNKNOWN_TABLE 0
#define ENERGY_TABLE 1
#define WAVEVECTOR_TABLE 2
#define WAVELENGTH_TABLE 3
#define FLUX_ADAPT_SET 0
#define FLUX_ADAPT_MULT 1
#define FLUX_ADAPT_ADD 2
FILE *hfile; /* id for data file */
char line[LINE_MAX_LENGTH];
long line_count_in_file = 0;
long line_count_in_array = 0;
long malloc_size = 100;
char flag_exit_loop = 0;
char flag_in_array = 0;
char Type_Table = UNKNOWN_TABLE;
double X,P;
double *Weight_Table;
double *KE_Table;
long Length_Table=0;
char Mode_Table = FLUX_ADAPT_SET;
char verbose=0;
int i;
double v2,K,L,E,new_p,slope;
double X1,X2,Y1,Y2;
/* end of declare */
%}
INITIALIZE
%{
KE_Table = (double*)malloc(malloc_size*sizeof(double));
Weight_Table = (double*)malloc(malloc_size*sizeof(double));
hfile = fopen(file, "r");
if(!hfile)
{
fprintf(stderr, "Error: %s : could not open input file '%s'\n", mccompcurname, file);
}
else /* now look for the data */
{
/* initialize data array */
if (strstr(options," k")
|| strstr(options," K ")
|| strstr(options,"wavevector"))
Type_Table = WAVEVECTOR_TABLE;
if (strstr(options,"omega")
|| strstr(options," e ")
|| strstr(options," E ")
|| strstr(options,"energy"))
Type_Table = ENERGY_TABLE;
if (strstr(options,"lambda")
|| strstr(options,"wavelength")
|| strstr(options," L "))
Type_Table = WAVELENGTH_TABLE;
if (strstr(options,"set"))
Mode_Table = FLUX_ADAPT_SET;
if (strstr(options,"add"))
Mode_Table = FLUX_ADAPT_ADD;
if (strstr(options,"multiply"))
Mode_Table = FLUX_ADAPT_MULT;
if (strstr(options,"verbose"))
verbose = 1;
/* do main loop */
while (!flag_exit_loop)
{
if (fgets(line, LINE_MAX_LENGTH, hfile) != NULL)
{ /* tries to read some informations */
line_count_in_file++;
for (i=0; (i < strlen(line)) && (i < LINE_MAX_LENGTH); i++) { line[i] = tolower(line[i]); }
if (strstr(line," k ")
|| strstr(line,"wavevector")) /* overrides options */
Type_Table = WAVEVECTOR_TABLE;
if (strstr(line," e ")
|| strstr(line,"omega")
|| strstr(line,"energy"))
Type_Table = ENERGY_TABLE;
if (strstr(line,"lambda")
|| strstr(line," l ")
|| strstr(line,"wavelength"))
Type_Table = WAVELENGTH_TABLE;
if (strstr(line,"set"))
Mode_Table = FLUX_ADAPT_SET;
if (strstr(line,"multiply"))
Mode_Table = FLUX_ADAPT_MULT;
if (strstr(line,"add"))
Mode_Table = FLUX_ADAPT_ADD;
if (strstr(line,"unactivate"))
Type_Table = UNKNOWN_TABLE;;
if (strstr(line,"verbose"))
verbose = 1;
/* tries to read 2 numbers */
if (sscanf(line,"%lg %lg",&X,&P) == 2)
{
/* if succeed and not in array : initialize array */
if (!flag_in_array)
{
flag_in_array = 1;
line_count_in_array = 0;
malloc_size = 0;
}
/* if succeed and in array : add (and realloc if necessary) */
if (line_count_in_array+100 > malloc_size)
{
malloc_size += 100;
KE_Table = (double*)realloc(KE_Table,malloc_size*sizeof(double));
Weight_Table = (double*)realloc(Weight_Table,malloc_size*sizeof(double));
}
KE_Table[line_count_in_array] = X;
Weight_Table[line_count_in_array] = P;
line_count_in_array++;
}
else
/* if does not succeed : set 'not in array' flag */
{
flag_in_array = 0;
}
}
else
flag_exit_loop = 1;
/* else : end of file */
}
Length_Table = line_count_in_array;
if (Length_Table < 2) Type_Table = UNKNOWN_TABLE; /* not enough points */
if (verbose)
{
printf("Flux : %i points in %s\n",Length_Table, file);
printf("Flux : data is [ ");
if (Type_Table == ENERGY_TABLE) printf("Energy");
if (Type_Table == WAVEVECTOR_TABLE) printf("Wavevector");
if (Type_Table == WAVELENGTH_TABLE) printf("Wavelength");
if (Type_Table == UNKNOWN_TABLE) printf("UNKNOWN (not used)");
printf(", Flux ]");
if (Mode_Table == FLUX_ADAPT_MULT) printf(" in multiply mode");
printf("\n");
}
fclose(hfile);
}
%}
TRACE
%{
PROP_Z0;
if (Type_Table && (x>xmin && x<xmax && y>ymin && y<ymax))
{
v2 = (vx*vx + vy*vy + vz*vz);
K = V2K*sqrt(v2); /* k */
L = (2*PI/K); /* lambda */
E = VS2E*v2; /* energy */
if (Type_Table == ENERGY_TABLE) X=E;
if (Type_Table == WAVEVECTOR_TABLE) X=K;
if (Type_Table == WAVELENGTH_TABLE) X=L;
/* table look up */
/* look for number just after X in table */
line_count_in_array=1;
while ((line_count_in_array <= Length_Table-1) && (KE_Table[line_count_in_array] < X)) line_count_in_array++;
X2 = KE_Table[line_count_in_array];
X1 = KE_Table[line_count_in_array-1];
Y2 = Weight_Table[line_count_in_array];
Y1 = Weight_Table[line_count_in_array-1];
if (X2-X1)
{
/* linear interpolation */
slope = (Y2-Y1)/(X2-X1);
new_p = Y1+slope*(X-X1);
if (Mode_Table == FLUX_ADAPT_MULT) p *= new_p;
else p = new_p;
}
else
ABSORB;
}
else
if (Type_Table) ABSORB;
%}
FINALLY
%{
free(KE_Table);
free(Weight_Table);
%}
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:
*
* none (see Source_Optimizer.comp)
*
*******************************************************************************/
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*3))
|| (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.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) passing in the Source_Optimizer, and then records a Reference source
* (step 2) as well as the state (at Source_Optimizer position) of neutrons
* reaching Monitor. The optimized source is defined as a fraction of the
* Reference source plus the distribution of 'good' neutrons reaching the
* Monitor. The optimization then starts (step 3), and focuses new neutrons on
* the Monitor_Optimizer. In fact it changes 'bad' neutrons into 'good' ones
* (that reach the Monitor). The overall flux is kept during process.
*
* 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
* ('absorb').
*
* Source_optimizer can be placed just after the source (for instance).
* 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). Overridden by 'auto' option.
* 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 sets 'step' as required.
* '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")
* A good optimization needs to record enough non optimized neutrons on Monitor
* during step 2. Typical enhancement in computation speed is by a factor 20.
*
*******************************************************************************/
/* 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 very efficient
Sep 29 1999 : v0.02 re-direct 'bad' neutrons instead of ABSORB (rand generator for nothing)
Oct 06 1999 : v0.03 installed options, corrected bugs, improved efficiency
Oct 21 1999 : v0.04 optim can be choosen for xy,v,s,p
*/
/* NOTE : The 'Optim_' variables are reserved by those components. */
/* other options : setxy, setv, setp, sets to precise what paremeters should be optimized */
/* default is 'setxy'+setv+sets' */
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 /* what was modified in last optim */
#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
#define OPTIM_DO_XY 1 /* what to optimize */
#define OPTIM_DO_V 2
#define OPTIM_DO_S 4
#define OPTIM_DO_P 8
/* 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 =1; /* 1 : continuous Source optimization */
char Optim_Flag_Recycle =0; /* record of neutron state changes by OPTIM_MOD_xx */
char Optim_Flag_KeepE =1; /* 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 */
char Optim_Flag_Type = OPTIM_DO_XY|OPTIM_DO_V|OPTIM_DO_S;
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,"set")) Optim_Flag_Type = 0;
if (strstr(options,"setxy")) Optim_Flag_Type |= OPTIM_DO_XY;
if (strstr(options,"setv")) Optim_Flag_Type |= OPTIM_DO_V;
if (strstr(options,"sets")) Optim_Flag_Type |= OPTIM_DO_S;
if (strstr(options,"setp")) Optim_Flag_Type |= OPTIM_DO_P;
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_Flag_Type & OPTIM_DO_V) && (Optim_Passing_vz[Optim_index] >= Optim_Source_vz[Optim_index]))
{ /* distribution achieved : redirect neutron near last neutron characteristic */
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_Flag_Type & OPTIM_DO_V) && (Optim_Passing_vx[Optim_index] >= Optim_Source_vx[Optim_index]))
{
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_Flag_Type & OPTIM_DO_V) && (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_Flag_Type & OPTIM_DO_XY) && (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_Flag_Type & OPTIM_DO_XY) && (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_Flag_Type & OPTIM_DO_S) && (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_Flag_Type & OPTIM_DO_S) && (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_Flag_Type & OPTIM_DO_P) && (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) */
if (Optim_Flag_Type & OPTIM_DO_V)
{
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;
}
if (Optim_Flag_Type & OPTIM_DO_XY)
{
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;
}
if (Optim_Flag_Type & OPTIM_DO_S)
{
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;
}
if (Optim_Flag_Type & OPTIM_DO_P)
{
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: ");
if (Optim_Flag_Continuous) fprintf(hfile,"continuous "); else fprintf(hfile,"fixed ");
if (Optim_Flag_Auto) fprintf(hfile,"auto ");
if (Optim_Flag_Absorb) fprintf(hfile,"absorb ");
if (Optim_Flag_KeepE) fprintf(hfile,"keepE "); else fprintf(hfile,"freeE ");
if (Optim_Flag_Verbose) fprintf(hfile,"verbose ");
if (Optim_Flag_Type & OPTIM_DO_XY) fprintf(hfile,"setxy ");
if (Optim_Flag_Type & OPTIM_DO_V) fprintf(hfile,"setv ");
if (Optim_Flag_Type & OPTIM_DO_S) fprintf(hfile,"sets ");
if (Optim_Flag_Type & OPTIM_DO_P) fprintf(hfile,"setp ");
fprintf(hfile,"\n");
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
More information about the mcstas-users
mailing list