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

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