Changeset 5110 for trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx

Timestamp:

09/17/04 18:36:25 (20 years ago)

Author:

moralejo

Message:

 Commited to CVS all changes to adapt headers to current C++ syntaxis

File:

• trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx (modified) (7 diffs)

trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx

@@ -4 +4 @@
 // (c) 2000 D. Petry 
 //
+// 15/09/2004, A. Moralejo: 
+// - Adapted to gcc 3.2 under root 3.05.07
+// - Fixed algorithm to calculate director cosines of incident photons.
+//   Former algorithm resulted in mirrored images on the camera (which we
+//   must not have after reflection on a parabollic dish).
+//
 /////////////////////////////////////////////////////////////////////////
 
@@ -9 +15 @@
 #include "starfield.h"
 
-#define PROGRAMID "$Id: starfield.cxx,v 1.3 2000-09-21 10:31:36 harald Exp $"
+#define PROGRAMID "$Id: starfield.cxx,v 1.4 2004-09-17 17:36:25 moralejo Exp $"
 
 /////////////////////////////////////////////////////////////////////////
 
-main(int argc, char **argv)
+int main(int argc, char **argv)
 {
 
@@ -32 +38 @@
   float lmin_nm[4] = {ULMIN_nm, BLMIN_nm, VLMIN_nm, RLMIN_nm}; // wave band definitions
   float lmax_nm[4] = {ULMAX_nm, BLMAX_nm, VLMAX_nm, RLMAX_nm}; 
-  float theta_rad, costheta, phi_rad, randtime, lambda_nm, xdum_m, ydum_m;
+  float theta_rad, costheta, sintheta, phi_rad, randtime, lambda_nm, xdum_m, ydum_m;
+  float cosa, cosA, sinA;
   float angdist;
 
@@ -88 +95 @@
   cout << "<Goddard Space Flight Center, Flight Dynamics Division (1998)>\n";
 
-  angdist = fmod( pars.catalog_fov_deg / cos( pars.ct_dec_rad ), 360.);
+  angdist = fmod( (float) (pars.catalog_fov_deg/cos( pars.ct_dec_rad )), 
+                  (float) 360.);
 
   if(angdist > 180.){ // too near to the pole, have to loop over all files
@@ -192 +200 @@
     }
     
-    stars[i].u = -1. * cos( stars[i].dec_rad ) * sin( stars[i].ra_rad - pars.ct_ra_rad ) / costheta; 
-
-    stars[i].v = -1. * ( sin( pars.ct_dec_rad ) * cos( stars[i].dec_rad ) * 
-                 cos( stars[i].ra_rad - pars.ct_ra_rad ) -
-                 cos( pars.ct_dec_rad ) * sin( stars[i].dec_rad )
-                 ) / costheta;
-    
+    sintheta = sqrt(1.-costheta*costheta);
+
+    //
+    // A. Moralejo, 15/09/2004  
+    // We want the director cosines of the down-going versors along the photon 
+    // incident directions. This is what the Reflector program expects (also from the 
+    // normal Corsika output). We have used simple spherical trigonometry to obtain
+    // the formulae.
+    //
+    
+    // "cosa" is the cosine of the angle "a" between the star direction and the direction
+    // defined by declination = 0  and right ascension = ct_ra. "a" is one of the sides of 
+    // a spherical triangle. It is a quantity needed for the calculations.
+
+    cosa = cos( stars[i].ra_rad - pars.ct_ra_rad ) * cos( stars[i].dec_rad );
+
+    // cosA is the angle between the great circle of constant right ascension = ct_ra  and
+    // the great circle defined by the star direction and the telescope direction. "A" would
+    // be the angle opposite to the side "a" of a spherical triangle (see above)
+
+    cosA = (cosa - costheta*cos(pars.ct_dec_rad)) / (sintheta*sin(pars.ct_dec_rad));
+
+    // We now want A to be defined in the 0 - 2pi range, so that it becomes the azimuth angle
+    // of the star in a system defined by the telescope direction. "A" will be defined between
+    // 0 and pi if sin (star_ra - ct_ra) > 0, and between pi and 2 pi otherwise.
+
+    sinA = sqrt (1.-cosA*cosA);
+    if ( sin(stars[i].ra_rad - pars.ct_ra_rad ) < 0. )
+      sinA *= -1.;
+
+    stars[i].u = -1. * sintheta * cosA;
+    stars[i].v = -1. * sintheta * sinA;
+    
+
+    // Old implementation, commented out 15/09/2004:
+    //
+    // This produced director cosines which, when fed to reflector, makes on the camera plane 
+    // a mirror-inverted image of the FOV (with respect to what one would see "by eye"). That 
+    // is NOT what we want! After reflection on the parabollic mirror, the image on the camera 
+    // is NOT mirrored, but just rotated by 180 degree. The new implementation above produces
+    // the correct FOV (tested with Reflector 0.6)
+    //
+    //    stars[i].u = -1. * cos( stars[i].dec_rad ) * sin( stars[i].ra_rad - pars.ct_ra_rad ) / costheta; 
+    //    stars[i].v = -1. * ( sin( pars.ct_dec_rad ) * cos( stars[i].dec_rad ) * 
+    //           cos( stars[i].ra_rad - pars.ct_ra_rad ) -
+    //           cos( pars.ct_dec_rad ) * sin( stars[i].dec_rad )
+    //           ) / costheta;
+    
+
     // calculate the "zenith angle" theta and "azimuth" phi of the star assuming 
     // the telecope points at the zenith 
@@ -211 +261 @@
       phi_rad = 2.*PI - acos(stars[i].u);
     }
-    
+
     // calculate number of photons
     
@@ -289 +339 @@
 
   convertcorsika( 
-                 ((int)pars.ct_ra_h*10)*1000 + (int)(abs(pars.ct_dec_deg)*10), // the file id
+                 ((int)pars.ct_ra_h*10)*1000 + (int)(fabs(pars.ct_dec_deg)*10), // the file id
                  totalphotinside, photons, pars.integtime_s, pars.verbose, pars.output_file);