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

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