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

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