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