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

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