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

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