1 gaskelld 1.1 subroutine trip_thru_target (narm, zpos, energy, theta, Eloss, radlen,
2 & mass, typeflag)
3
4 implicit none
5 include 'simulate.inc'
6
7 integer narm
8 integer typeflag !1=generate eloss, 2=min, 3=max, 4=most probable
9 real*8 zpos, energy, mass, theta
10 real*8 Eloss, radlen
11 real*8 forward_path, side_path
12 real*8 s_target, s_Al, s_kevlar, s_air, s_mylar ! distances travelled
13 real*8 Eloss_target, Eloss_Al,Eloss_air ! energy losses
14 real*8 Eloss_kevlar,Eloss_mylar ! (temporary)
15 real*8 z_can,t,atmp,btmp,ctmp,costmp,th_can !for the pudding-can target.
16 logical liquid
17
18 s_Al = 0.0
19 liquid = targ.Z.lt.2.4
20
21 if (abs(zpos) .gt. (targ.length/2.+1.d-5)) then
22 gaskelld 1.1 write(6,*) 'call to trip_thru_target has |zpos| > targ.length/2.'
23 write(6,*) 'could be numerical error, or could be error in target offset'
24 write(6,*) 'zpos=',zpos,' targ.length/2.=',targ.length/2.
25 endif
26
27 ! Which particle are we interested in?
28
29 goto (10,20,30) narm
30
31 ! The incoming electron
32
33 10 continue
34 s_target = (targ.length/2. + zpos) / abs(cos(targ.angle))
35 if (liquid) then !liquid target
36 if (targ.can .eq. 1) then !beer can (2.8 mil endcap)
37 s_Al = s_Al + 0.0028*inch_cm
38 else if (targ.can .eq. 2) then !pudding can (5 mil Al, for now)
39 s_Al = s_Al + 0.0050*inch_cm
40 endif
41 endif
42
43 gaskelld 1.1 ! ... compute distance in radiation lengths and energy loss
44 radlen = s_target/targ.X0_cm + s_Al/X0_cm_Al
45 call enerloss_new(s_target,targ.rho,targ.Z,targ.A,energy,mass,
46 & typeflag,Eloss_target)
47 call enerloss_new(s_Al,rho_Al,Z_Al,A_Al,energy,mass,typeflag,Eloss_Al)
48 Eloss = Eloss_target + Eloss_Al
49 return
50
51 ! The scattered electron
52
53 ! ... compute distances travelled assuming HMS (SOS)
54 ! ...... 16.0 (8.0) mil of Al-foil for the target chamber exit foil
55 ! ......~15 cm air between scattering chamber and spectrometer vacuum
56 ! ...... 15.0 (5.0) mil of kevlar + 5.0 (3.0) mil mylar for the
57 ! ...... spectrometer entrance foil
58 ! ...... ASSUMES liquid targets are 2.65" wide, and have 5.0 mil Al side walls.
59
60 ! ... compute distances travelled assuming HRS (E. Schulte thesis)
61 ! ...... 13.0 mil of Al-foil for the target chamber exit foil
62 ! ...... WILD GUESS: ~15 cm air between scattering chamber and HRS vacuum
63 ! ...... 10.0 mil of Kapton for spectrometer entrance (Use mylar, since
64 gaskelld 1.1 ! ..... X0=28.6cm for Kapton, X0=28.7cm for Mylar)
65 ! ...... ASSUMES liquid targets are 2.65" wide, and have 5.0 mil Al side walls.
66
67 ! ... For SHMS, use SOS windows for now (just a space filler for now).
68
69 20 continue
70 if (electron_arm.eq.1) then !electron is in HMS
71 s_Al = 0.016*inch_cm
72 s_air = 15
73 s_kevlar = 0.015*inch_cm
74 s_mylar = 0.005*inch_cm
75 forward_path = (targ.length/2.-zpos) / abs(cos(theta+targ.angle))
76 else if (electron_arm.eq.2) then !SOS
77 s_Al = 0.008*inch_cm
78 s_air = 15
79 s_kevlar = 0.005*inch_cm
80 s_mylar = 0.003*inch_cm
81 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
82 else if (electron_arm.eq.3) then !HRS-R
83 s_Al = 0.013*inch_cm
84 s_air = 15
85 gaskelld 1.1 s_kevlar = 0.*inch_cm
86 s_mylar = 0.010*inch_cm
87 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
88 else if (electron_arm.eq.4) then !HRS-L
89 s_Al = 0.013*inch_cm
90 s_air = 15
91 s_kevlar = 0.*inch_cm
92 s_mylar = 0.010*inch_cm
93 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
94 else if (electron_arm.eq.5 .or. electron_arm.eq.6) then !SHMS
95 s_Al = 0.008*inch_cm
96 s_air = 15
97 s_kevlar = 0.005*inch_cm
98 s_mylar = 0.003*inch_cm
99 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
100 endif
101 s_target = forward_path
102
103 if (liquid) then
104 if (targ.can .eq. 1) then !beer can
105 side_path = 1.325*inch_cm / abs(sin(theta))
106 gaskelld 1.1 if (forward_path.lt.side_path) then
107 s_Al = s_Al + 0.005*inch_cm / abs(cos(theta))
108 else
109 s_target = side_path
110 s_Al = s_Al + 0.005*inch_cm / abs(sin(theta))
111 endif
112 else if (targ.can .eq. 2) then !pudding can (5 mil Al, for now)
113
114 ! this is ugly. Solve for z position where particle intersects can. The
115 ! pathlength is then (z_intersect - z_scatter)/cos(theta)
116 ! Angle from center to z_intersect is acos(z_intersect/R). Therefore the
117 ! angle between the particle and wall is pi/2 - (theta - theta_intersect)
118
119 t=tan(theta)**2
120 atmp=1+t
121 btmp=-2*zpos*t
122 ctmp=zpos**2*t-(targ.length/2.)**2
123 z_can=(-btmp+sqrt(btmp**2-4.*atmp*ctmp))/2./atmp
124 side_path = (z_can - zpos)/abs(cos(theta))
125 s_target = side_path
126 costmp=z_can/(targ.length/2.)
127 gaskelld 1.1 if (abs(costmp).le.1) then
128 th_can=acos(z_can/(targ.length/2.))
129 else if (abs(costmp-1.).le.0.000001) then
130 th_can=0. !extreme_trip_thru_target can give z/R SLIGHTLY>1.0
131 else
132 stop 'z_can > can radius in target.f !!!'
133 endif
134 s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
135 endif
136
137 endif
138
139 ! ... compute distance in radiation lengths and energy loss
140 radlen = s_target/targ.X0_cm + s_Al/X0_cm_Al + s_air/X0_cm_air +
141 > s_kevlar/X0_cm_kevlar + s_mylar/X0_cm_mylar
142 call enerloss_new(s_target,targ.rho,targ.Z,targ.A,energy,mass,
143 & typeflag,Eloss_target)
144 call enerloss_new(s_Al,rho_Al,Z_Al,A_Al,energy,mass,typeflag,Eloss_Al)
145 call enerloss_new(s_air,rho_air,Z_air,A_air,energy,mass,typeflag,
146 & Eloss_air)
147 call enerloss_new(s_kevlar,rho_kevlar,Z_kevlar,A_kevlar,energy,
148 gaskelld 1.1 & mass,typeflag,Eloss_kevlar)
149 call enerloss_new(s_mylar,rho_mylar,Z_mylar,A_mylar,energy,mass,
150 & typeflag,Eloss_mylar)
151 Eloss = Eloss_target + Eloss_Al + Eloss_air + Eloss_kevlar + Eloss_mylar
152
153 return
154
155 ! The scattered proton
156
157 30 continue
158 if (hadron_arm.eq.1) then !proton in HMS
159 s_Al = 0.016*inch_cm
160 s_air = 15
161 s_kevlar = 0.015*inch_cm
162 s_mylar = 0.005*inch_cm
163 forward_path = (targ.length/2.-zpos) / abs(cos(theta+targ.angle))
164 else if (hadron_arm.eq.2) then !SOS
165 s_Al = 0.008*inch_cm
166 s_air = 15
167 s_kevlar = 0.005*inch_cm
168 s_mylar = 0.003*inch_cm
169 gaskelld 1.1 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
170 else if (hadron_arm.eq.3) then !HRS-R
171 s_Al = 0.013*inch_cm
172 s_air = 15
173 s_kevlar = 0.*inch_cm
174 s_mylar = 0.010*inch_cm
175 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
176 else if (hadron_arm.eq.4) then !HRS-L
177
178 ! if (abs(zpos/targ.length).lt.0.00001) write(6,*) 'SOMEONE LEFT THE SAFETY WINDOW IN SIMC!!!!!!!'
179 ! s_Al = 0.125*inch_cm
180
181 s_Al = 0.013*inch_cm
182 s_air = 15
183 s_kevlar = 0.*inch_cm
184 s_mylar = 0.010*inch_cm
185 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
186 else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then !SHMS
187 s_Al = 0.008*inch_cm
188 s_air = 15
189 s_kevlar = 0.005*inch_cm
190 gaskelld 1.1 s_mylar = 0.003*inch_cm
191 forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
192 endif
193
194 s_target = forward_path
195 if (liquid) then
196 if (targ.can .eq. 1) then !beer can
197 side_path = 1.325*inch_cm/ abs(sin(theta))
198 if (forward_path.lt.side_path) then
199 s_Al = s_Al + 0.005*inch_cm / abs(cos(theta))
200 else
201 s_target = side_path
202 s_Al = s_Al + 0.005*inch_cm / abs(sin(theta))
203 endif
204 else if (targ.can .eq. 2) then !pudding can (5 mil Al, for now)
205
206 ! this is ugly. Solve for z position where particle intersects can. The
207 ! pathlength is then (z_intersect - z_scatter)/cos(theta)
208 ! Angle from center to z_intersect is acos(z_intersect/R). Therefore the
209 ! angle between the particle and wall is pi/2 - (theta - theta_intersect)
210
211 gaskelld 1.1 t=tan(theta)**2
212 atmp=1+t
213 btmp=-2*zpos*t
214 ctmp=zpos**2*t-(targ.length/2.)**2
215 z_can=(-btmp+sqrt(btmp**2-4.*atmp*ctmp))/2./atmp
216 side_path = (z_can - zpos)/abs(cos(theta))
217 s_target = side_path
218 costmp=z_can/(targ.length/2.)
219 if (abs(costmp).le.1) then
220 th_can=acos(z_can/(targ.length/2.))
221 else if (abs(costmp-1.).le.0.000001) then
222 th_can=0. !extreme_trip_thru_target can give z/R SLIGHTLY>1.0
223 else
224 stop 'z_can > can radius in target.f !!!'
225 endif
226 s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
227 endif
228
229 endif
230
231 ! ... compute energy losses
232 gaskelld 1.1
233 radlen = s_target/targ.X0_cm + s_Al/X0_cm_Al + s_air/X0_cm_air +
234 > s_kevlar/X0_cm_kevlar + s_mylar/X0_cm_mylar
235 call enerloss_new(s_target,targ.rho,targ.Z,targ.A,energy,mass,
236 & typeflag,Eloss_target)
237 call enerloss_new(s_Al,rho_Al,Z_Al,A_Al,energy,mass,typeflag,Eloss_Al)
238 call enerloss_new(s_air,rho_air,Z_air,A_air,energy,mass,typeflag,
239 & Eloss_air)
240 call enerloss_new(s_kevlar,rho_kevlar,Z_kevlar,A_kevlar,energy,mass,
241 & typeflag,Eloss_kevlar)
242 call enerloss_new(s_mylar,rho_mylar,Z_mylar,A_mylar,energy,mass,
243 & typeflag,Eloss_mylar)
244
245 Eloss=Eloss_target+Eloss_Al+Eloss_air+Eloss_kevlar+Eloss_mylar
246
247 return
248 end
249
250 !-----------------------------------------------------------------
251
252 subroutine extreme_trip_thru_target(ebeam, the, thp, pe, pp, z, m)
253 gaskelld 1.1
254 implicit none
255 include 'target.inc'
256 include 'constants.inc'
257 structure /limits/
258 real*8 min, max
259 end structure
260
261 record /limits/ the, thp, pe, pp, z, betap
262 integer i
263 real*8 th_corner, th_corner_min, th_corner_max
264 real*8 E1, E2, E3, E4, t1, t2, t3, t4
265 real*8 zz, th1, th2, m
266 real*8 ebeam, energymin, energymax
267 logical liquid
268
269 real*8 zero
270 parameter (zero=0.0d0) !double precision zero for subroutine calls
271
272 !Given limiting values for the electron/proton angles, the z-position in the
273 !target, and beta for the proton, determine min and max losses in target (and
274 gaskelld 1.1 !corresponding amounts of material traversed) This procedure requires knowledge
275 !of the shape of the target, of course, that's why I've put it in here!
276 !Bothering to compute this at all, given that it's going to be used for slop
277 !limits, just proves what a glowing example of Obsessive Compulsive Disorder at
278 !work i truly am.
279
280 liquid = targ.Z.lt.2.4
281
282 ! Incoming electron
283 call trip_thru_target(1, z.max, ebeam, zero, targ.Eloss(1).max,
284 & targ.teff(1).max, Me, 3)
285 call trip_thru_target(1, z.min, ebeam, zero, targ.Eloss(1).min,
286 & targ.teff(1).min, Me, 2)
287
288 ! Scattered electron
289 C Above a few MeV, the energy loss increases as a function of electron
290 C energy, so for any energies we're dealing the max(min) energy loss should
291 C correspond to the max(min) spectrometer energy. Note that this isn't
292 C the perfect range, but it's easier than reproducing the generated limits here
293
294 energymin=pe.min
295 gaskelld 1.1 energymax=pe.max
296 ! ... solid target is infinitely wide
297 if (.not.liquid) then
298 call trip_thru_target(2, z.min, energymax, the.max,
299 & targ.Eloss(2).max, targ.teff(2).max, Me, 3)
300 call trip_thru_target(2, z.max, energymin, the.min,
301 & targ.Eloss(2).min, targ.teff(2).min, Me, 2)
302
303 ! ... liquid
304 else if (targ.can .eq. 1) then !beer can
305 if (z.max.ge.targ.length/2.) then
306 th_corner_max = target_pi/2.
307 else
308 th_corner_max = atan(1.25*inch_cm/(targ.length/2.-z.max))
309 endif
310 th_corner_min = atan(1.25*inch_cm/(targ.length/2.-z.min))
311
312 ! ... max loss: do we have access to a corner shot? try a hair on either
313 ! ... side, front or side walls could be thicker (too lazy to check!)
314 if (th_corner_min.le.the.max .and. th_corner_max.ge.the.min) then
315 th_corner = max(th_corner_min,the.min)
316 gaskelld 1.1 zz = targ.length/2. - 1.25*inch_cm/tan(th_corner)
317 th1 = th_corner-.0001
318 th2 = th_corner+.0001
319 else
320 zz = z.min
321 th1 = the.min
322 th2 = the.max
323 endif
324 call trip_thru_target(2, zz, energymax, th1, E1, t1, Me, 3)
325 call trip_thru_target(2, zz, energymax, th2, E2, t2, Me, 3)
326 targ.Eloss(2).max = max(E1,E2)
327 targ.teff(2).max = max(t1,t2)
328
329 ! ........ min loss: try all possibilities (in min case, no way
330 ! ........ an intermediate z or th will do the trick).
331 targ.Eloss(2).min = 1.d10
332 do i = 0, 3
333 call trip_thru_target (2, z.min+int(i/2)*(z.max-z.min), energymin,
334 > the.min+mod(i,2)*(the.max-the.min), E1, t1, Me, 2)
335 if (E1 .lt. targ.Eloss(2).min) then
336 targ.Eloss(2).min = E1
337 gaskelld 1.1 targ.teff(2).min = t1
338 endif
339 enddo
340
341 else if (targ.can .eq. 2) then !pudding can
342
343 ! ... for pudding can, max loss occurs at lowest scattering angle, where
344 ! ... the outgoing particle leaves the can at z=0. Therefore, take minimum
345 ! ... scattering angle, project back from z=0, x=radius to get z_init.
346 ! ... Minimum z_init is -radius.
347 zz = -(targ.length/2.)/tan(the.min)
348 zz = max (zz,(-targ.length/2.))
349 call trip_thru_target(2, zz, energymax, the.min, targ.Eloss(2).max,
350 & targ.teff(2).max, Me, 3)
351
352 ! ........ min loss: try all possibilities (in min case, no way
353 ! ........ an intermediate z or th will do the trick). Actually, this
354 ! ........ may not be true for the pudding can target, but I'm leaving the
355 ! ........ code alone, for now.
356 targ.Eloss(2).min = 1.d10
357 do i = 0, 3
358 gaskelld 1.1 call trip_thru_target (2, z.min+int(i/2)*(z.max-z.min), energymin,
359 > the.min+mod(i,2)*(the.max-the.min), E1, t1, Me, 2)
360 if (E1 .lt. targ.Eloss(2).min) then
361 targ.Eloss(2).min = E1
362 targ.teff(2).min = t1
363 endif
364 enddo
365 endif
366
367 ! Scattered proton. As you can see I'm sufficiently lazy to make
368 ! the code work out whether high or low beta
369 ! maximizes or minimizes loss ... always lower beta --> higher
370 ! losses, it seems
371
372 ! ... compute extreme energy values
373 energymin = sqrt(pp.min**2 + m**2)
374 energymax = sqrt(pp.max**2 + m**2)
375 ! ... compute extreme betap values
376 betap.min = pp.min / sqrt(pp.min**2 + m**2)
377 betap.max = pp.max / sqrt(pp.max**2 + m**2)
378 ! ... solid target is infinitely wide
379 gaskelld 1.1 if (.not.liquid) then
380 call trip_thru_target (3, z.min, energymin, thp.max, E1, t1, m, 3)
381 call trip_thru_target (3, z.min, energymax, thp.max, E2, t2, m, 3)
382 targ.Eloss(3).max = max(E1,E2)
383 targ.teff(3).max = max(t1,t2)
384
385 call trip_thru_target (3, z.max, energymin, thp.min, E1, t1, m, 2)
386 call trip_thru_target (3, z.max, energymax, thp.min, E2, t2, m, 2)
387 targ.Eloss(3).min = min(E1,E2)
388 targ.teff(3).min = min(t1,t2)
389
390 ! ... liquid
391 else if (targ.can .eq. 1) then !beer can
392 if (z.max .ge. targ.length/2.) then
393 th_corner_max = target_pi/2.
394 else
395 th_corner_max = atan(1.25*inch_cm/(targ.length/2.-z.max))
396 endif
397 th_corner_min = atan(1.25*inch_cm/(targ.length/2.-z.min))
398 ! ........ max loss: do we have access to a corner shot?
399 ! ........ try a hair on either side, front or side walls could be
400 gaskelld 1.1 ! thicker (too lazy to check!)
401 if (th_corner_min.le.thp.max .and. th_corner_max.ge.thp.min) then
402 th_corner = max(th_corner_min,thp.min)
403 zz = targ.length/2. - 1.25*inch_cm/tan(th_corner)
404 th1 = th_corner-.0001
405 th2 = th_corner+.0001
406 else
407 zz = z.min
408 th1 = thp.min
409 th2 = thp.max
410 endif
411 call trip_thru_target (3, zz, energymin, th1, E1, t1, m, 3)
412 call trip_thru_target (3, zz, energymin, th2, E2, t2, m, 3)
413 call trip_thru_target (3, zz, energymax, th1, E3, t3, m, 3)
414 call trip_thru_target (3, zz, energymax, th2, E4, t4, m, 3)
415 targ.Eloss(3).max = max(E1,E2,E3,E4)
416 targ.teff(3).max = max(t1,t2,t3,t4)
417
418 ! ........ min loss: try all possibilities (in min case, no way
419 ! ........ an intermediate z or th will do the trick)
420 targ.Eloss(3).min = 1.d10
421 gaskelld 1.1 do i = 0, 3
422 call trip_thru_target (3, z.min+int(i/2)*(z.max-z.min), energymin,
423 > thp.min+mod(i,2)*(thp.max-thp.min), E1, t1, m, 2)
424 if (E1 .lt. targ.Eloss(3).min) then
425 targ.Eloss(3).min = E1
426 targ.teff(3).min = t1
427 zz = z.min+int(i/2)*(z.max-z.min)
428 th1 = thp.min+mod(i,2)*(thp.max-thp.min)
429 endif
430 enddo
431 call trip_thru_target (3, zz, energymax, th1, E1, t1, m, 2)
432 targ.Eloss(3).min = min(targ.Eloss(3).min, E1)
433 else if (targ.can .eq. 2) then !pudding can
434
435 ! ... for pudding can, max loss occurs at lowest scattering angle, where
436 ! ... the outgoing particle leaves the can at z=0. Therefore, take minimum
437 ! ... scattering angle, project back from z=0, x=radius to get z_init.
438 ! ... Minimum z_init is -radius.
439 zz = -(targ.length/2.)/tan(the.min)
440 zz = max (zz,(-targ.length/2.))
441
442 gaskelld 1.1 call trip_thru_target (3, zz, energymin, thp.min, E1, t1, m, 3)
443 call trip_thru_target (3, zz, energymax, thp.min, E2, t2, m, 3)
444 targ.Eloss(3).max = max(E1,E2)
445 targ.teff(3).max = max(t1,t2)
446
447 ! ........ min loss: try all possibilities (in min case, no way
448 ! ........ an intermediate z or th will do the trick). This may be
449 ! ........ wrong for the pudding can targ. Check it out later.
450 targ.Eloss(3).min = 1.d10
451 do i = 0, 3
452 call trip_thru_target (3, z.min+int(i/2)*(z.max-z.min), energymin,
453 > thp.min+mod(i,2)*(thp.max-thp.min), E1, t1, m, 2)
454 if (E1 .lt. targ.Eloss(3).min) then
455 targ.Eloss(3).min = E1
456 targ.teff(3).min = t1
457 zz = z.min+int(i/2)*(z.max-z.min)
458 th1 = thp.min+mod(i,2)*(thp.max-thp.min)
459 endif
460 enddo
461 call trip_thru_target (3, zz, energymax, th1, E1, t1, m, 2)
462 targ.Eloss(3).min = min(targ.Eloss(3).min, E1)
463 gaskelld 1.1
464 endif
465
466 *JRA*! Extreme multiple scattering
467 *JRA*
468 *JRA*! ........ don't consider multiple scattering of incoming electron
469 *JRA* targ.musc_max(1) = 0.
470
471 ! Extreme multiple scattering. Use nominal beam energy rather than minimum
472 ! (should be close enough)
473
474 call extreme_target_musc(ebeam,1.d0,
475 > targ.teff(1).max,targ.musc_max(1),targ.musc_nsig_max)
476 call extreme_target_musc(pe.min,1.d0,
477 > targ.teff(2).max,targ.musc_max(2),targ.musc_nsig_max)
478 call extreme_target_musc(pp.min,betap.min,
479 > targ.teff(3).max,targ.musc_max(3),targ.musc_nsig_max)
480
481 return
482 end
483
484 gaskelld 1.1 !------------------------------------------------------------------
485
486 subroutine target_musc(p,beta,teff,dangles)
487
488 implicit none
489
490 real*8 Es, epsilon, nsig_max
491 parameter (Es = 13.6) !MeV
492 parameter (epsilon = 0.088)
493 parameter (nsig_max = 3.5)
494
495 real*8 p, beta, teff, dangles(2), dangle, r
496 real*8 theta_sigma
497 real*8 gauss1
498
499 if (p.lt.25.) write(6,*)
500 > 'Momentum passed to target_musc should be in MeV, but p=',p
501
502 ! Compute rms value for planar scattering angle distribution, cf. PDB
503 ! Note teff is thickness of material, in radiation lengths.
504
505 gaskelld 1.1 theta_sigma = Es/p/beta * sqrt(teff) * (1+epsilon*log10(teff))
506
507 ! Compute scattering angles in perpendicular planes.
508 ! Generate two Gaussian numbers BELOW nsig_max.
509
510 dangles(1) = theta_sigma * gauss1(nsig_max)
511 dangles(2) = theta_sigma * gauss1(nsig_max)
512
513 return
514
515 ! Return info about extreme multiple scattering
516
517 entry extreme_target_musc(p, beta, teff, dangle, r)
518
519 theta_sigma = Es/p/beta * sqrt(teff) * (1+epsilon*log10(teff))
520 dangle = theta_sigma * nsig_max
521 r = nsig_max
522 return
523 end
|