source: trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c@ 1614

Last change on this file since 1614 was 1614, checked in by bigongia, 22 years ago
*** empty log message ***
File size: 22.6 KB
Line 
1#include <stdio.h>
2#include <math.h>
3#include "diag.h"
4#include "init.h"
5#include "lagrange.h"
6/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
7
8/* random numbers */
9#define RandomNumber ranf()
10
11/* Speed of Light in vacuum, in m/s */
12#define Speed_of_Light_vacuum 299792458.0f
13#define Speed_of_Light_air (Speed_of_Light_vacuum / 1.000293f)
14
15/* Speed of Light in vacuum, in cm/ns */
16#define Speed_of_Light_vacuum_cmns (Speed_of_Light_vacuum / 1.0e7f)
17#define Speed_of_Light_air_cmns (Speed_of_Light_air / 1.0e7f)
18
19/* Macros */
20#define SQR(A) ((A)*(A))
21#define NORM(A) ((float) sqrt((SQR(A[0]))+(SQR(A[1]))+(SQR(A[2]))))
22
23/* Function declarations */
24extern float ranf(void);
25void rnormal(double *r, int n);
26
27/* Static definitions */
28float OmegaCT[3][3];
29float OmegaICT[3][3];
30float Omega[3][3];
31float OmegaI[3][3];
32
33static double NormalRandomNumbers[500];
34/*
35 From photons on ground, i.e. observation level,
36 to photons on focal plane, i.e. chamber !
37
38 Mirror reflectivity
39 Mirror reflection
40 Photon position on chamber
41 Position smearing
42 Timing
43
44 Returned values:
45
46 0 OK photon reached the chamber
47
48 1 Photon lost due to mirror reflectivity
49 2 Photon lost because out of mirror
50 3 Photon lost due to black spot
51 4 Photon lost because reflected out of chamber
52
53*/
54
55int ph2cph(photon *ph, cphoton *cph)
56{
57
58 float u, v, w; /* photon director cosines */
59
60 float r[3]; /* photon trajectory */
61 float x[3]; /* position of the photon on ground */
62
63 float rCT[3]; /* photon trajectory in the system of the CT */
64 float xCT[3]; /* photon position on ground (CT) */
65 float rm[3]; /* photon trajectory in the system of a mirror */
66 float xmm[3]; /* intermediate values */
67 float xm[3]; /* photon position on ground */
68 float xcut[3]; /* location of the cut sphere-trajectory */
69 float xcutCT[3]; /* location of the cut sphere-trajectory (CT) */
70
71 float rnor[3], rnorm; /* normal in that point */
72
73 float rrefl[3]; /* reflected vector, from the mirror */
74 float rreflCT[3]; /* reflected vector, from the CT */
75 float xcam[3]; /* where the photon hits the camera plane */
76
77 float calpha; /* cos(alpha=angle incident/normal) */
78 float phi; /* angle between photon and camera plane */
79
80 float a, b, c, t; /* intermediate variables */
81
82 float d; /* minimum distance trajectory-mirror center */
83
84 float wl; /* photon wavelength */
85 float reflec; /* reflectivity for a photon */
86
87 float h; /* photon production height */
88
89 int i, k; /* simple counters */
90 int i_mirror=-1; /* number of a given mirror */
91
92
93 float distmirr, distmirr2; /* distances used in MAGIC reflection routine */
94 float sx, sy;
95 float t1, t2;
96 float dummy = 0.;
97
98
99 void makeOmega(float theta, float phi);
100 void makeOmegaI(float theta, float phi);
101 void applyMxV(float M[3][3], float *V, float *Vp);
102 float Lin2Curv(float x);
103
104 /* begin code */
105
106 /* get photon wawelength */
107
108 wl = ph->w;
109
110 /* get position on ground */
111
112 x[0] = ph->x;
113 x[1] = ph->y;
114 x[2] = 0.0; /* ground => obs. level => z=0 */
115
116 /* get director cosines x,y on ground */
117
118 r[0] = ph->u;
119 r[1] = ph->v;
120
121// AM 11/2002: fixed line below: u v are the direction cosines of the
122// *downgoing* photon. Hence, third component must be negative!
123// This was a serious bug affecting all versions before 0.6 (see TDAS
124// note on Reflector program 0.6).
125
126 r[2] = (float) -sqrt(1.0 - r[0]*r[0] - r[1]*r[1]);
127
128 /* get photon time and production height */
129
130 h = ph->h;
131
132 /* CBC */
133 Debug("@0 x r %f %f %f %f %f %f\n", x[0], x[1], x[2],
134 r[0], r[1], r[2]);
135
136 /*
137 x[0] = 125.0;
138 x[1] = 125.0;
139 x[2] = 0.0; */ /* ground => obs. level => z=0 */
140
141 /* get director cosines x,y on ground */
142
143 /*
144 r[0] = 0.0;
145 r[1] = 0.0;
146 r[2] = 1.0;
147 */
148 /* CBC */
149
150 /*!@'
151
152 @#### Reflectivity of the mirrors.
153
154 We make a 3rd. order interpolation using Lagrange
155 polynomials, in order to calculate the reflectivity of the
156 mirror for that wavelength. Further developments will
157 include also a incidence-angle dependence (which is not very
158 important).
159
160 */
161
162 /* ++
163 FILTER: REFLECTIVITY R(lambda)
164 -- */
165
166 /* find data point to be used in Lagrange interpolation (-> k) */
167
168 FindLagrange(Reflectivity,k,wl);
169
170 /* if random > reflectivity then goes to the TOP of the loop again */
171
172 reflec = Lagrange(Reflectivity,k,wl);
173
174 if ( RandomNumber > reflec ) return 1;
175
176
177 /*!@'
178
179 @#### Reflection on mirrors.
180 We calculate reflected photon direction
181
182 */
183
184 /* ++
185 REFLECTION
186 -- */
187
188 Debug("@1 x r %f %f %f %f %f %f\n", x[0], x[1], x[2], r[0], r[1], r[2]);
189
190 /* change to the system of the CT */
191
192 applyMxV( OmegaCT, x, xCT );
193 applyMxV( OmegaCT, r, rCT );
194
195 /* CBC */
196 Debug("@2 xCT rCT %f %f %f %f %f %f\n", xCT[0], xCT[1], xCT[2],
197 rCT[0], rCT[1], rCT[2]);
198 /* CBC */
199
200
201 /*
202 before moving to the system of the mirror
203 we look whether the photon hits a mirror or not
204
205 calculate the intersection of the trajectory of the photon
206 with the GLOBAL DISH !!!
207 we reproduce the calculation of the coefficients of the
208 second order polynomial in z (=xCT[2]), made with
209 Mathematica
210 */
211
212 /*
213 * In[1]:= parab:=z-(x^2+y^2)/(4F)
214 * par1=parab /. {x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
215 *
216 * Out[1]=
217 * u (z - z0) 2 v (z - z0) 2
218 * (x0 + ----------) + (y0 + ----------)
219 * w w
220 * z - ---------------------------------------
221 * 4 F
222 *
223 * In[2]:= CoefficientList[ExpandAll[par1*4F*w^2],z]
224 *
225 * Out[2]=
226 * 2 2 2 2
227 * {-(w x0 ) - w y0 + 2 u w x0 z0 +
228 *
229 * 2 2 2 2
230 * 2 v w y0 z0 - u z0 - v z0 ,
231 *
232 * 2 2
233 * 4 F w - 2 u w x0 - 2 v w y0 + 2 u z0 +
234 *
235 * 2 2 2
236 * 2 v z0, -u - v }
237 */
238
239 /* the z coordinate is calculated */
240
241 a = - SQR(rCT[0]) - SQR(rCT[1]);
242
243 b = (float) (4.0*ct_Focal_mean*SQR(rCT[2])
244 - 2.0*rCT[0]*rCT[2]*xCT[0] - 2.0*rCT[1]*rCT[2]*xCT[1]
245 + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2]);
246
247 /* FIXED Lines below, May 2002, AM : formerly (up to V0.4)
248 * there was a confusion between telescope coordinates xCT and
249 * the original coordinates x. Thanks to T. Hengstebeck for
250 * reporting the bug.
251 */
252
253 c = 2*rCT[0]*rCT[2]*xCT[0]*xCT[2] + 2*rCT[1]*rCT[2]*xCT[1]*xCT[2]
254 - SQR(rCT[2])*SQR(xCT[0]) - SQR(rCT[2])*SQR(xCT[1])
255 - SQR(rCT[0])*SQR(xCT[2]) - SQR(rCT[1])*SQR(xCT[2]);
256
257
258 if ( fabs(a) < 1.e-6 ) {
259
260 /* only one value */
261
262 xcut[2] = -c / b;
263
264 } else {
265
266 /* Introduce positiveness check, AM 3/7/2002 */
267
268 dummy = b*b - 4.0*a*c;
269
270 if (dummy < 0.) /* No intersection */
271 return 2;
272
273 d = (float) sqrt(dummy);
274
275 /* two possible values for z */
276
277 t1 = (float) ((-b+d) / (2.0*a));
278 t2 = (float) ((-b-d) / (2.0*a));
279
280 /* z must be the minimum of t1 and t2 */
281
282 xcut[2] = (t1 < t2) ? t1 : t2;
283
284 }
285
286 /*
287 xcut[] is NOW the cut between the GLOBAL dish of MAGIC and
288 the trajectory of the photon
289 */
290
291 xcut[0] = xCT[0] + rCT[0]/rCT[2]*(xcut[2]-xCT[2]);
292 xcut[1] = xCT[1] + rCT[1]/rCT[2]*(xcut[2]-xCT[2]);
293
294
295 /* CBC */
296 Debug("@3 xcut %f %f\n", xcut[0], xcut[1]);
297 /* CBC */
298
299 /* convert to Curvilinear distance over the parabolic dish */
300
301 sx = Lin2Curv( xcut[0] );
302 sy = Lin2Curv( xcut[1] );
303
304 /* CBC */
305 Debug("@4 sx sy %f %f\n", sx, sy);
306 /* CBC */
307
308 /* is it outside the dish? */
309
310 if ((fabs(sx) > 850.0) || (fabs(sy) > 850.0)) {
311 /*
312 cout << "CONDITION 1 !" << endl << flush;
313 cout << '1';
314 */
315 return 2;
316 }
317
318 /* calculate the mirror to be used */
319
320 distmirr = 1000000.0f;
321
322 for (i=0; i<ct_NMirrors && distmirr>=ct_RMirror; ++i) {
323 distmirr2 = (float) sqrt(SQR(ct_data[i].x - xcut[0]) +
324 SQR(ct_data[i].y - xcut[1]) +
325 SQR(ct_data[i].z - xcut[2]));
326 if (distmirr2 < distmirr) {
327 i_mirror = i;
328 distmirr = distmirr2;
329 }
330 }
331
332 /*
333 the mirror to use is i_mirror (calculated several lines above)
334 check whether the photon is outside the nearest (this) mirror
335 */
336
337 if ((fabs(ct_data[i_mirror].sx - sx) > ct_RMirror) ||
338 (fabs(ct_data[i_mirror].sy - sy) > ct_RMirror)) {
339 /*
340 cout << "CONDITION 2 !" << endl << flush;
341 cout << '2';
342 */
343 return 2;
344 }
345
346 /* CBC */
347 Debug("@5 theta phi %f %f\n", ct_data[i_mirror].theta,
348 ct_data[i_mirror].phi);
349 /* CBC */
350
351 /* calculate matrices for the mirror */
352
353 makeOmega (-ct_data[i_mirror].theta,
354 ct_data[i_mirror].phi);
355 makeOmegaI(-ct_data[i_mirror].theta,
356 ct_data[i_mirror].phi);
357
358 /* change to the system of the mirror */
359
360 /* CBC */
361 Debug("@6 mirror %f %f %f\n",ct_data[i_mirror].x,
362 ct_data[i_mirror].y,
363 ct_data[i_mirror].z);
364 /* CBC */
365
366 /* first translation... */
367
368 xmm[0] = xCT[0] - ct_data[i_mirror].x;
369 xmm[1] = xCT[1] - ct_data[i_mirror].y;
370 xmm[2] = xCT[2] - ct_data[i_mirror].z;
371
372
373 /* CBC */
374 Debug("@7 xmm %f %f %f\n", xmm[0], xmm[1], xmm[2]);
375 /* CBC */
376
377 /* ...then rotation */
378
379 applyMxV( Omega, xmm, xm );
380 applyMxV( Omega, rCT, rm );
381
382 /* CBC */
383 Debug("@8 xm rm %f %f %f %f %f %f\n", xm[0], xm[1], xm[2],
384 rm[0], rm[1], rm[2]);
385 /* CBC */
386
387 /*
388 the vector rCT should be normalized, and
389 so the vector rm remains normalized as well, but, anyhow...
390 */
391
392 rnorm = NORM( rm );
393 rm[0] /= rnorm;
394 rm[1] /= rnorm;
395 rm[2] /= rnorm;
396
397 /* CBC */
398 Debug("@9 rm-norm %f %f %f\n", rm[0], rm[1], rm[2]);
399 /* CBC */
400
401 /*
402 calculate the intersection of the trajectory of the photon
403 with the mirror
404 we reproduce the calculation of the coefficients of the
405 second order polynomial in z (=xm[2]), made with
406 Mathematica
407 */
408
409 /*
410 * In[1]:= esfera:=x^2+y^2+(z-R)^2-R^2;
411 * recta:={x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
412 *
413 * In[2]:= esfera
414 *
415 * 2 2 2 2
416 * Out[2]= -R + x + y + (-R + z)
417 *
418 * In[3]:= recta
419 *
420 * u (z - z0) v (z - z0)
421 * Out[3]= {x -> x0 + ----------, y -> y0 + ----------}
422 * w w
423 *
424 * In[4]:= esf=esfera /. recta
425 *
426 * 2 2 u (z - z0) 2 v (z - z0) 2
427 * Out[4]= -R + (-R + z) + (x0 + ----------) + (y0 + ----------)
428 * w w
429 *
430 * In[5]:= coefs=CoefficientList[ExpandAll[esf],z]
431 *
432 * 2 2 2 2
433 * 2 2 2 u x0 z0 2 v y0 z0 u z0 v z0
434 * Out[5]= {x0 + y0 - --------- - --------- + ------ + ------,
435 * w w 2 2
436 * w w
437 *
438 * 2 2 2 2
439 * 2 u x0 2 v y0 2 u z0 2 v z0 u v
440 * > -2 R + ------ + ------ - ------- - -------, 1 + -- + --}
441 * w w 2 2 2 2
442 * w w w w
443 * In[6]:= Simplify[ExpandAll[coefs*w^2]]
444 *
445 * 2 2 2 2 2 2
446 * Out[6]= {w (x0 + y0 ) - 2 w (u x0 + v y0) z0 + (u + v ) z0 ,
447 *
448 * 2 2 2 2 2
449 * > -2 (R w - u w x0 + u z0 + v (-(w y0) + v z0)), u + v + w }
450 *
451 */
452
453 /*
454 the z coordinate is calculated, using the coefficient
455 shown above
456 */
457
458 a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]);
459
460 b = (float) (-2*(2.*ct_data[i_mirror].f*SQR(rm[2])
461 - rm[0]*rm[2]*xm[0]
462 + SQR(rm[0])*xm[2]
463 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2])));
464
465 c = (SQR(rm[2])*(SQR(xm[0]) + SQR(xm[1]))
466 - 2*rm[2]*(rm[0]*xm[0] + rm[1]*xm[1])*xm[2]
467 + (SQR(rm[0]) + SQR(rm[1]))*SQR(xm[2]));
468
469 d = (float) sqrt( b*b - 4.0*a*c );
470
471 /* two possible values for z */
472
473 t1 = (float) ((-b+d) / (2.0*a));
474 t2 = (float) ((-b-d) / (2.0*a));
475
476 /* z must be the minimum of t1 and t2 */
477
478 xcut[2] = (t1 < t2) ? t1 : t2;
479 xcut[0] = xm[0] + rm[0]/rm[2]*(xcut[2]-xm[2]);
480 xcut[1] = xm[1] + rm[1]/rm[2]*(xcut[2]-xm[2]);
481
482 /* CBC */
483 Debug("@10 xcut %f %f %f\n", xcut[0], xcut[1], xcut[2]);
484 /* CBC */
485
486 /*
487 ++
488 BLACK SPOTS: If the photon hits the black spot, it's lost
489 --
490 */
491
492 if ( sqrt(SQR(xcut[0]) + SQR(xcut[1])) < ct_BlackSpot_rad ) {
493 /*
494 cout << "CONDITION 3!\n" << flush;
495 cout << '3';
496 */
497 return 3;
498 }
499
500 /*
501 if we still have the photon, we continue with the reflexion;
502 we calculate normal vector in this point and normalize:
503 */
504
505 rnor[0] = 2.0f*xcut[0];
506 rnor[1] = 2.0f*xcut[1];
507 rnor[2] = (float) (2.0*(xcut[2] - 2.0*ct_data[i_mirror].f));
508
509 /* CBC */
510 Debug("@11 rnor %f %f %f\n", rnor[0], rnor[1], rnor[2]);
511 /* CBC */
512
513// Changed AM, 11/2002: now we use the normal vector going "outwards"
514// from inside the sphere (=removed minus sign in normalization below).
515// It is easier to do so, since now the vector rm indicating the
516// photon direction also goes from the front to the back of the mirror.
517
518 rnorm = NORM( rnor );
519 rnor[0] /= rnorm;
520 rnor[1] /= rnorm;
521 rnor[2] /= rnorm;
522
523 /* CBC */
524 Debug("@12 rnor-norm %f %f %f\n", rnor[0], rnor[1], rnor[2]);
525 /* CBC */
526
527 /*
528 now, both "normal" vector and original trajectory are
529 normalized
530 just project the original vector in the normal, and
531 take it as the "mean" position of the original and
532 the "reflected" vector
533 from this, we can calculate the "reflected" vector
534 calpha = cos(angle(rnor,rm))
535 */
536
537// AM 11/2002: removed absolute value in scalar
538// product below (it is now unnecessary):
539
540 calpha = (float) (rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
541
542 /* CBC */
543 Debug("@13 calpha %f\n", calpha);
544 /* CBC */
545
546 /* finally!!! we have the reflected trajectory of the photon */
547
548
549 rrefl[0] = (float) (2.0*rnor[0]*calpha - rm[0]);
550 rrefl[1] = (float) (2.0*rnor[1]*calpha - rm[1]);
551 rrefl[2] = (float) (2.0*rnor[2]*calpha - rm[2]);
552
553 /* CBC */
554 Debug("@14 rrefl %f %f %f\n", rrefl[0], rrefl[1], rrefl[2]);
555 /* CBC */
556
557 rnorm = NORM( rrefl );
558 rrefl[0] /= rnorm;
559 rrefl[1] /= rnorm;
560 rrefl[2] /= rnorm;
561
562 /* CBC */
563 Debug("@15 rrefl-norm %f %f %f\n", rrefl[0], rrefl[1], rrefl[2]);
564 /* CBC */
565
566 /* let's go back to the coordinate system of the CT */
567
568 /* first rotation... */
569
570 applyMxV( OmegaI, xcut, xcutCT);
571 applyMxV( OmegaI, rrefl, rreflCT);
572
573 /* CBC */
574 Debug("@16 xcutCT rreflCT %f %f %f %f %f %f\n", xcutCT[0], xcutCT[1],
575 xcutCT[2], rreflCT[0], rreflCT[1], rreflCT[2]);
576 /* CBC */
577
578 /* ...then translation */
579
580 xcutCT[0] += ct_data[i_mirror].x;
581 xcutCT[1] += ct_data[i_mirror].y;
582 xcutCT[2] += ct_data[i_mirror].z;
583
584 /* CBC */
585 Debug("@17 xcutCT %f %f %f\n", xcutCT[0], xcutCT[1], xcutCT[2]);
586 /* CBC */
587
588 /*
589 calculate intersection of this trajectory and the camera plane
590 in the system of the CT, this plane is z = ct_Focal
591 */
592
593 t = (ct_Focal_mean - xcutCT[2]) / rreflCT[2];
594
595 xcam[0] = xcutCT[0] + rreflCT[0]*t;
596 xcam[1] = xcutCT[1] + rreflCT[1]*t;
597 xcam[2] = xcutCT[2] + rreflCT[2]*t;
598
599 /* CBC */
600 Debug("@18 xcam %f %f %f\n", xcam[0], xcam[1], xcam[2]);
601 /* CBC */
602
603 /*
604 ++
605 AXIS DEVIATION: We introduce it here just as a first order
606 correction, by modifying the position of the reflected photon.
607 --
608 */
609
610 xcam[0] += AxisDeviation[0][i_mirror];
611 xcam[1] += AxisDeviation[1][i_mirror];
612
613 /* CBC */
614 Debug("@19 xcam-AD %f %f \n", xcam[0], xcam[1]);
615 /* CBC */
616
617 /*
618 ++
619 SMEARING: We apply the point spread function for the mirrors
620 --
621 */
622
623 /* get two N(0;1) random numbers */
624
625 rnormal( NormalRandomNumbers, 2 );
626
627 /* modify the Cphoton position in the camera */
628
629 xcam[0] += (float) (NormalRandomNumbers[0] * ct_PSpread_mean);
630 xcam[1] += (float) (NormalRandomNumbers[1] * ct_PSpread_mean);
631
632 /* CBC */
633 Debug("@20 xcam-SM %f %f \n", xcam[0], xcam[1]);
634 /* CBC */
635
636 /* check whether the photon goes out of the camera */
637
638 if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) {
639 return 4;
640 }
641
642 /*
643 ++
644 ANGLE OF INCIDENCE
645 --
646
647 calculate angle of incidence between tray. and camera plane
648 the camera plane is
649 0 x + 0 y + z - ct_Focal_mean = 0 => (A,B,C,D) = (0,0,1,-ct_Focal_mean)
650 from Table 3.20 "Tasch. der Math."
651 */
652
653 phi = (float) asin(rreflCT[2]);
654
655 /*
656 ++
657 TIMING
658 --
659 */
660
661 /* calculate the new time of the photon (in the camera) */
662
663 t = ph->t;
664
665 /*
666 substract path from the mirror till the ground, 'cos
667 the photon actually hit the mirror!!
668 */
669
670 t = (float) (t + ((( xm[2] > 0. ) ? -1.0 : +1.0) *
671 sqrt( SQR(xm[0] - xcut[0]) +
672 SQR(xm[1] - xcut[1]) +
673 SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns));
674
675 /* add path from the mirror till the camera */
676
677 t = (float) (t + sqrt( SQR(xcutCT[0] - xcam[0]) +
678 SQR(xcutCT[1] - xcam[1]) +
679 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns);
680
681 /* show it */
682
683 Debug("@22 %f %f %f\n"
684 "@23 %f %f %f %f %f %f\n"
685 "@24 %f %f %d %f %f %f %f\n"
686 "@25 %f %f %f %f\n\n",
687 xCT[0], xCT[1], xCT[2], rCT[0], rCT[1], rCT[2],
688 xcut[0], xcut[1], xcut[2],
689 sx, sy, i_mirror, ct_data[i_mirror].sx, ct_data[i_mirror].sy,
690 ct_data[i_mirror].sx - sx, ct_data[i_mirror].sy - sy,
691 xcam[0], xcam[1], xcam[2], phi);
692
693 /* Output */
694
695 /* AM Nov 2002: Added one further change of coordinates so that the camera
696 * images have the "right" orientation: they will appear as seen by an
697 * observer on ground, standing behind the mirror dish and looking towards
698 * the camera. Formerly cph->x and cph->y were simply xcam[0] and xcam[1].
699 */
700
701 cph->x = -xcam[1];
702 cph->y = -xcam[0];
703
704 cph->u = r[0];
705 cph->v = r[1];
706 cph->t = t;
707 cph->h = h;
708 cph->phi = phi;
709
710 return 0;
711
712} /* end of ph2cph */
713
714
715/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
716
717!---------------------------------------------------------------------
718 @name makeOmega
719
720 @desc function to calculate the matrix Omega(theta,phi)
721
722 @var theta Angle theta of the transformation
723 @var phi Angle phi of the transformation
724
725 @date Sat Jun 27 05:58:56 MET DST 1998
726----------------------------------------------------------------------
727 @function
728*/
729
730void
731makeOmega (float theta, float phi)
732{
733 static float ct, st, cp, sp;
734
735 /* shortcuts for cosine and sine of theta and phi */
736 ct = (float) cos(theta);
737 st = (float) sin(theta);
738 cp = (float) cos(phi);
739 sp = (float) sin(phi);
740
741 /* save values in the array (see top of file) */
742 Omega[0][0] = cp*ct;
743 Omega[0][1] = sp*ct;
744 Omega[0][2] = -st;
745
746 Omega[1][0] = -sp;
747 Omega[1][1] = cp;
748 Omega[1][2] = 0;
749
750 Omega[2][0] = cp*st;
751 Omega[2][1] = sp*st;
752 Omega[2][2] = ct;
753}
754
755
756/*
757!---------------------------------------------------------------------
758 @name makeOmegaI
759
760 @desc function to calculate the matrix Omega-1(theta,phi)
761
762 @var theta Angle theta of the transformation
763 @var phi Angle phi of the transformation
764
765 @date Sat Jun 27 05:58:56 MET DST 1998
766----------------------------------------------------------------------
767 @function
768*/
769
770void
771makeOmegaI(float theta, float phi)
772{
773 static float ct, st, cp, sp;
774
775 /* shortcuts for cosine and sine of theta and phi */
776 ct = (float) cos(theta);
777 st = (float) sin(theta);
778 cp = (float) cos(phi);
779 sp = (float) sin(phi);
780
781 /* save values in the array (see top of file) */
782 OmegaI[0][0] = cp*ct;
783 OmegaI[0][1] = -sp;
784 OmegaI[0][2] = cp*st;
785
786 OmegaI[1][0] = sp*ct;
787 OmegaI[1][1] = cp;
788 OmegaI[1][2] = sp*st;
789
790 OmegaI[2][0] = -st;
791 OmegaI[2][1] = 0;
792 OmegaI[2][2] = ct;
793}
794
795
796/*
797!---------------------------------------------------------------------
798 @name applyMxv
799
800 @desc returns the vector v' such that v' = M x v
801
802 @var M matrix of the transformation
803 @var v vector to be multiplied
804 @var vi resulting vector
805
806 @date Sat Jun 27 05:58:56 MET DST 1998
807----------------------------------------------------------------------
808 @function
809*/
810
811void
812applyMxV(float M[3][3], float *V, float *Vp)
813{
814 Vp[0] = (M[0][0] * V[0] +
815 M[0][1] * V[1] +
816 M[0][2] * V[2]);
817 Vp[1] = (M[1][0] * V[0] +
818 M[1][1] * V[1] +
819 M[1][2] * V[2]);
820 Vp[2] = (M[2][0] * V[0] +
821 M[2][1] * V[1] +
822 M[2][2] * V[2]);
823}
824
825/*
826!---------------------------------------------------------------------
827 @name Lin2Curv
828
829 @desc Linear (Euclidean) to Curvilinear distance
830
831 @var x Radial distance from the axis of the paraboloid
832
833 @return Curvilinear distance over the parabolic shape
834
835 @date Wed Jul 8 15:25:39 MET DST 1998
836----------------------------------------------------------------------
837 @function
838*/
839
840float
841Lin2Curv(float x)
842{
843 x /= 100.f;
844 return ((x + 0.000144175317185f * x * x * x)*100.f);
845}
846
847/*!---------------------------------------------------------------------
848// @name rnormal
849//
850// @desc returns n(=2k) normaly distributed numbers
851//
852// @var *r pointer to a vector where we write the numbers
853// @var n how many numbers do we generate
854//
855// @date Sat Jun 27 05:58:56 MET DST 1998
856//----------------------------------------------------------------------
857// @function */
858
859void rnormal(double *r, int n)
860{
861
862 double z1, z2;
863 int i;
864
865 for (i=0; i<n; i+=2) {
866
867 z1 = RandomNumber;
868 z2 = RandomNumber;
869
870 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
871 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
872
873 }
874
875}
Note: See TracBrowser for help on using the repository browser.