1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
/*******************************************************************************
*
* Component: Conics_EH
*
* %Identification
* Written by: Peter Wilendrup and Erik Knudsen <br>(derived from Giacomo Resta skeleton-component)
* Origin: DTU
* Release: McStas 2.7
* Date: September 2021
*
* %Description
* Implements a set of nshells Wolter Ellipsoid/Hyperboloid pairs using conics.h from ConicTracer.
*
* The component has two distinct modes of specifying the geometry:
* a) Via the radii vector, parametrized from largest to smallest radius with a length of 'nshells'
*
* b) By specifying radii rmax and rmin, between which a quadratic law distributes 'nshells' surfaces.
*
* The mirrors are assumed to be touching at the mid-optic plane, i.e. there is no gap between primary
* and secondary mirror.
* By definition the ratio between primary and secondary mirror glancing angles is 1/3.
* At present a single m-value is used for all mirrors.
*
* %Parameters
*
* Input parameters:
* nshells: [1] Number of nested shells to expect
* rmin: [m] Midoptic plane radius of innermost mirror pair.
* rmax: [m] Midoptic plane radius of outermost mirror pair.
* radii: [m] Optional vector of radii (length should match nshells)
* m: [1] Critical angle of mirrors as multiples of Ni_c.
* focal_length_u: [m] Focal length (upstream) of the mirror pairs.
* focal_length_d: [m] Focal length (downstream) of the mirror pairs.
* le: [m] Paraboloid mirror length.
* lh: [m] Hyperboloid mirror length.
* disk: [ ] Flag. If nonzero, insert a disk to block the central area within the innermost mirror.
* Qc: [AA-1] Critical scattering vector
* R0: [1] Reflectivity at Low angles for reflectivity curve approximation
* alpha: [AA] Slope of reflectivity for reflectivity curve approximation
* W: [AA-1] Width of supermirror cut-off
* transmit: [1] Fraction of statistics to assign to transmitted beam - NOT YET IMPLEMENTED
* mirr_thick: [m] Thickness of mirror shell surfaces - NOT YET IMPLEMENTED
* %End
*
*******************************************************************************/
DEFINE COMPONENT Conics_EH
SETTING PARAMETERS (
rmin=0.0031416, rmax=0.05236,focal_length_u=10, focal_length_d=10, le=0.25, lh=0.25,
int nshells=4, m=1, mirr_thick=0, int disk=1, vector radii=NULL,
R0 = 0.99, Qc = 0.021, W = 0.003, alpha = 6.07, transmit = 0)
SHARE
%{
%include "ref-lib"
%include "conic.h"
%include "read_table-lib"
%}
DECLARE
%{
//Scene where all geometry is added to
Scene s;
%}
INITIALIZE
%{
ConicSurf *pm;
double th_c, alpha_p, alpha_h, fp2, dr,rr, cH, theta_1, theta_2, theta_i;
int i;
s=makeScene();
/* Mode a, vector of radii */
if (radii) {
for (i=0;i<nshells;i++){
rr=radii[i];
Point pi = makePoint(0,rr,0);
pm=addEllipsoid(focal_length_d, -focal_length_u, pi, -le, 0, m, R0, Qc, W, alpha, &s);
addHyperboloid( focal_length_d, focal_length_d, pi, 0, lh, m, R0, Qc, W, alpha, &s);
}
} else {
/* Mode b, use quadratic law to distribute the shells */
dr=nshells>1?(rmax-rmin)/(nshells-1):0;
double constant, quadratic;
quadratic = (rmax-rmin)/(rmax*rmax - rmin*rmin);
constant = rmax - quadratic*rmax*rmax;
for (i=0;i<nshells;i++){
rr = rmax-dr*i;
rr = constant + quadratic*rr*rr; // Quadratic distribution of radius covers angular space better
//printf("--------------------------------------------------------------------");
printf("rr = %lf\n", rr);
//printf("--------------------------------------------------------------------");
Point pi = makePoint(0,rr,0);
//pm=addEllipsoid(-focal_length_u, focal_length_d , pi, -le, 0, m,R0,Qc,W,alpha,&s);
//addHyperboloid( focal_length_d, focal_length_d*2, pi, 0, lh, m,R0,Qc,W,alpha,&s);
theta_1 = atan(rr/focal_length_u);
theta_2 = atan(rr/focal_length_d);
theta_i = 0.25*(theta_1 + theta_2);
cH = fabs(0.5*(rr/tan(theta_2 - 2.0*theta_i) - focal_length_d));
pm=addEllipsoid(focal_length_d + 2.0*cH, -focal_length_u, pi, -le, 0, m, R0, Qc, W, alpha, &s);
addHyperboloid( focal_length_d, focal_length_d + 2.0*cH, pi, 0, lh, m, R0, Qc, W, alpha, &s);
}
}
if (disk) {
addDisk(pm->zs,0.0,rConic(pm->ze,*pm),&s);
}
%}
TRACE
%{
/* "_mctmp_a" defines a "silicon" state variable in underlying conic.h functions */
_mctmp_a=0;
traceSingleNeutron(_particle,s);
if (!_particle->_absorbed) {
SCATTER;
}
%}
FINALLY %{
//Mainly Writes Inline Detector Data
finishSimulation(&s);
%}
MCDISPLAY
%{
double zz = 0;
//Enlarge xy-plane when mcdisplay is ran with --zoom
magnify("xy");
//Draw xy-axis contour for Conic Surfaces
int i;
for (i = 0; i < s.num_c; i++) {
double step = (s.c[i].ze-s.c[i].zs)/100;
double cz;
int draw=-1;
for (cz = s.c[i].zs+step; cz <= s.c[i].ze; cz+= step) {
draw++;
double rp = rConic(cz-step,s.c[i]);
double rc = rConic(cz, s.c[i]);
double rx,ry;
int j;
double theta;
for (j = 0; j < 12; j++) {
theta = 2.0*PI*j/12.0;
rx = rp*cos(theta);
ry = rp*sin(theta);
line(rx,ry,cz-step-zz,rx,ry,cz-zz);
}
if (draw==0) {
circle("xy", 0, 0, cz-step-zz, rp);
}
if (draw==19) draw=-1;
}
}
//Draw xy-axis cross hairs for Disks
//Local variables to control maximal display-size of cross-hairs
for (i = 0; i < s.num_di; i++) {
double r0=s.di[i].r0;
double r1=s.di[i].r1;
double z0=s.di[i].z0;
if (r0>1.0) r0=1.0;
if (r1>1.0) r1=1.0;
line(r0, 0, z0-zz, r1, 0, z0-zz);
line(-r0, 0, z0-zz, -r1, 0, z0-zz);
line(0, r0, z0-zz, 0, r1,z0-zz);
line(0, -r0, z0-zz, 0, -r1,z0-zz);
}
%}
END
|
