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

Last change on this file since 1621 was 1621, checked in by bigongia, 22 years ago
*** empty log message ***
File size: 22.8 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 /* AM, 15/11/2002: changed sign to get the angle between photon trajectory
654 * and camera plane positive! This had to be changed because now the vector
655 * indicating the reflected photon direction has the opposite sign!
656 */
657
658 phi = (float) -asin(rreflCT[2]);
659
660 /*
661 ++
662 TIMING
663 --
664 */
665
666 /* calculate the new time of the photon (in the camera) */
667
668 t = ph->t;
669
670 /*
671 substract path from the mirror till the ground, 'cos
672 the photon actually hit the mirror!!
673 */
674
675 t = (float) (t + ((( xm[2] > 0. ) ? -1.0 : +1.0) *
676 sqrt( SQR(xm[0] - xcut[0]) +
677 SQR(xm[1] - xcut[1]) +
678 SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns));
679
680 /* add path from the mirror till the camera */
681
682 t = (float) (t + sqrt( SQR(xcutCT[0] - xcam[0]) +
683 SQR(xcutCT[1] - xcam[1]) +
684 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns);
685
686 /* show it */
687
688 Debug("@22 %f %f %f\n"
689 "@23 %f %f %f %f %f %f\n"
690 "@24 %f %f %d %f %f %f %f\n"
691 "@25 %f %f %f %f\n\n",
692 xCT[0], xCT[1], xCT[2], rCT[0], rCT[1], rCT[2],
693 xcut[0], xcut[1], xcut[2],
694 sx, sy, i_mirror, ct_data[i_mirror].sx, ct_data[i_mirror].sy,
695 ct_data[i_mirror].sx - sx, ct_data[i_mirror].sy - sy,
696 xcam[0], xcam[1], xcam[2], phi);
697
698 /* Output */
699
700 /* AM Nov 2002: Added one further change of coordinates so that the camera
701 * images have the "right" orientation: they will appear as seen by an
702 * observer on ground, standing behind the mirror dish and looking towards
703 * the camera. Formerly cph->x and cph->y were simply xcam[0] and xcam[1].
704 */
705
706 cph->x = -xcam[1];
707 cph->y = -xcam[0];
708
709 cph->u = r[0];
710 cph->v = r[1];
711 cph->t = t;
712 cph->h = h;
713 cph->phi = phi;
714
715 return 0;
716
717} /* end of ph2cph */
718
719
720/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
721
722!---------------------------------------------------------------------
723 @name makeOmega
724
725 @desc function to calculate the matrix Omega(theta,phi)
726
727 @var theta Angle theta of the transformation
728 @var phi Angle phi of the transformation
729
730 @date Sat Jun 27 05:58:56 MET DST 1998
731----------------------------------------------------------------------
732 @function
733*/
734
735void
736makeOmega (float theta, float phi)
737{
738 static float ct, st, cp, sp;
739
740 /* shortcuts for cosine and sine of theta and phi */
741 ct = (float) cos(theta);
742 st = (float) sin(theta);
743 cp = (float) cos(phi);
744 sp = (float) sin(phi);
745
746 /* save values in the array (see top of file) */
747 Omega[0][0] = cp*ct;
748 Omega[0][1] = sp*ct;
749 Omega[0][2] = -st;
750
751 Omega[1][0] = -sp;
752 Omega[1][1] = cp;
753 Omega[1][2] = 0;
754
755 Omega[2][0] = cp*st;
756 Omega[2][1] = sp*st;
757 Omega[2][2] = ct;
758}
759
760
761/*
762!---------------------------------------------------------------------
763 @name makeOmegaI
764
765 @desc function to calculate the matrix Omega-1(theta,phi)
766
767 @var theta Angle theta of the transformation
768 @var phi Angle phi of the transformation
769
770 @date Sat Jun 27 05:58:56 MET DST 1998
771----------------------------------------------------------------------
772 @function
773*/
774
775void
776makeOmegaI(float theta, float phi)
777{
778 static float ct, st, cp, sp;
779
780 /* shortcuts for cosine and sine of theta and phi */
781 ct = (float) cos(theta);
782 st = (float) sin(theta);
783 cp = (float) cos(phi);
784 sp = (float) sin(phi);
785
786 /* save values in the array (see top of file) */
787 OmegaI[0][0] = cp*ct;
788 OmegaI[0][1] = -sp;
789 OmegaI[0][2] = cp*st;
790
791 OmegaI[1][0] = sp*ct;
792 OmegaI[1][1] = cp;
793 OmegaI[1][2] = sp*st;
794
795 OmegaI[2][0] = -st;
796 OmegaI[2][1] = 0;
797 OmegaI[2][2] = ct;
798}
799
800
801/*
802!---------------------------------------------------------------------
803 @name applyMxv
804
805 @desc returns the vector v' such that v' = M x v
806
807 @var M matrix of the transformation
808 @var v vector to be multiplied
809 @var vi resulting vector
810
811 @date Sat Jun 27 05:58:56 MET DST 1998
812----------------------------------------------------------------------
813 @function
814*/
815
816void
817applyMxV(float M[3][3], float *V, float *Vp)
818{
819 Vp[0] = (M[0][0] * V[0] +
820 M[0][1] * V[1] +
821 M[0][2] * V[2]);
822 Vp[1] = (M[1][0] * V[0] +
823 M[1][1] * V[1] +
824 M[1][2] * V[2]);
825 Vp[2] = (M[2][0] * V[0] +
826 M[2][1] * V[1] +
827 M[2][2] * V[2]);
828}
829
830/*
831!---------------------------------------------------------------------
832 @name Lin2Curv
833
834 @desc Linear (Euclidean) to Curvilinear distance
835
836 @var x Radial distance from the axis of the paraboloid
837
838 @return Curvilinear distance over the parabolic shape
839
840 @date Wed Jul 8 15:25:39 MET DST 1998
841----------------------------------------------------------------------
842 @function
843*/
844
845float
846Lin2Curv(float x)
847{
848 x /= 100.f;
849 return ((x + 0.000144175317185f * x * x * x)*100.f);
850}
851
852/*!---------------------------------------------------------------------
853// @name rnormal
854//
855// @desc returns n(=2k) normaly distributed numbers
856//
857// @var *r pointer to a vector where we write the numbers
858// @var n how many numbers do we generate
859//
860// @date Sat Jun 27 05:58:56 MET DST 1998
861//----------------------------------------------------------------------
862// @function */
863
864void rnormal(double *r, int n)
865{
866
867 double z1, z2;
868 int i;
869
870 for (i=0; i<n; i+=2) {
871
872 z1 = RandomNumber;
873 z2 = RandomNumber;
874
875 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
876 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
877
878 }
879
880}
Note: See TracBrowser for help on using the repository browser.