source: trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx@ 18066

/////////////////////////////////////////////////////////////////////////
// Starfield Generator
//
// (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).
//
// 04/10/2004, A. Moralejo:
// - Added to Makefile  root-config --libs  to the libraries. Otherwise
//   the program does not link in some machines in Munich! 
//
/////////////////////////////////////////////////////////////////////////


#include "starfield.h"

#define PROGRAMID "$Id: starfield.cxx,v 1.5 2004-10-04 11:34:49 moralejo Exp $"

/////////////////////////////////////////////////////////////////////////

int main(int argc, char **argv)
{

  parameters pars;
  star stars[iMAXSTARS];
  photon *photons;  

  ifstream in; 
  FILE *catfile;
  char parfilename[160];
  char catfilename[160];

  
  int i, j, k, numstars, subnumstars, starnumber, totalnumphot, photinside;
  int totalphotinside, idum;
  int istart_ra_h, iend_ra_h;
  int nph[4]; // numbers of photons in the four wavebands
  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, sintheta, phi_rad, randtime, lambda_nm, xdum_m, ydum_m;
  float cosa, cosA, sinA;
  float angdist;

  //welcome

  cout << "This is STARFIELD. (c) 2000 D. Petry\n";
  cout << PROGRAMID << "\n";

  // check command line arguments
  
  if(argc == 1){
    sprintf(parfilename, "starfield.par");
  }
  else{ // a filename was given
       sprintf(parfilename, "%s", argv[1]);
  }

  // read parameter file

  in.open(parfilename, ios::in); // open the parameter file

  if(!in){
    cout << "Failed to open " << parfilename << "\n" 
         << "Value of stream \"in\" was " << in << "\n";
    cout << "\nThere shoud be a parameter file "<< parfilename
         <<" in the working directory with the following format:\n-----\n";
    pars.usage(&cout);
    cout << "-----\nExiting.\n";
    exit(1);
  }

  cout << "Opened " << parfilename << " for reading ...\n";
  
  if( !(pars.readparameters(&in)) ){ // read not OK?
    if(!in.eof()){ // was the error not due to EOF?
      cout << "Error: rdstate = " << in.rdstate() << "\n";
    }
    else{
      cout << "Error: premature EOF in parameter file.\n";
    }
    exit(1);
  }
  
  in.close();

  // Allocate memory for photons

  photons=new photon[iMAXPHOT];


  // prepare loop over star catalog files

  cout << "SKY2000 - Master Star Catalog - Star Catalog Database, Version 2\n";
  cout << "Sande C.B., Warren W.H.Jr., Tracewell D.A., Home A.T., Miller A.C.\n";
  cout << "<Goddard Space Flight Center, Flight Dynamics Division (1998)>\n";

  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

     istart_ra_h = 0;
     iend_ra_h = 23;

  }
  else{ // can loop over selected files only

    istart_ra_h = (int) (pars.ct_ra_h - angdist / 360. * 24.) - 1; 
    iend_ra_h = (int) (pars.ct_ra_h + angdist / 360. * 24. ) + 1; 

  }
  
  //read catalog

  i = 0;
  for (j = istart_ra_h; j <= iend_ra_h; j++){

    subnumstars = 0;

    if ( j <  0){
      idum = j + 24;
    }
    else if ( j > 23 ){
      idum = j - 24;
    }
    else {
        idum = j;
    }

    sprintf(catfilename, "%s/sky%02d.dat", pars.datapath, idum);

    if((catfile = fopen(catfilename, "r")) == NULL){ // open the star catalog
      cout << "Failed to open " << catfilename << "\n";
      exit(1);
    }
    cout << "Opened file " << catfilename << " for reading ...\n";
    
    while( stars[i].readstar(catfile, pars.verbose) ){ // read next star OK
      
      angdist = acos( cos( pars.ct_dec_rad ) * cos( stars[i].dec_rad ) * 
        cos( stars[i].ra_rad - pars.ct_ra_rad ) +
        sin( pars.ct_dec_rad ) * sin( stars[i].dec_rad ) );
      
      if( angdist < pars.catalog_fov_deg / 180. * PI ){ // star in Field Of View?

        stars[i].calcmissingmags(pars.verbose);
        if (pars.verbose) stars[i].printstar();
        
        if( stars[i].umag > -100. ){
          i++; // accept star
          subnumstars++;
        }
        else{
          cout << "Star rejected.\n";
        }
      
        if( i > iMAXSTARS ){
          i--;
          cout << "Error: Star memory full. Accepted " << i << " stars.\n";
          break;
        }
      }    
    }
    if( feof(catfile) ){ // was  EOF reached?
      cout << "EOF reached; accepted "<< subnumstars << " stars from this segment.\n";
    }
    if( ferror(catfile) ){ // did an error occur?
      cout << "Error while reading catalog file.\n";
      exit(1);
    }
    fclose(catfile);
    if(i == iMAXSTARS){
      break;
    }
  }

  cout << "Accepted "<< i << " stars in total.\n";
  numstars = i;

  // loop over all photons from all stars, filling their fields
  
  totalnumphot=0;

  totalphotinside=0; 
    
  for(i=0; i<numstars;i++){
    
    starnumber=stars[i].icatnum;
    
    // calculate director cosines (see Montenbruck & Pfleger, 1989, p. 196)
    
    costheta = cos( pars.ct_dec_rad ) * cos( stars[i].dec_rad ) * 
      cos( stars[i].ra_rad - pars.ct_ra_rad ) +
      sin( pars.ct_dec_rad ) * sin( stars[i].dec_rad );
    
    if(costheta == 0.){
      cout << "Star number " << i << " (catalog number " << stars[i].icatnum <<
        ") seems to be at 90 degrees distance from optical axis.\n ... will ignore it.";
      continue;
    }
    
    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 
    // take into account the ambiguity of acos() 
    
    theta_rad = acos(costheta);
    
    if( stars[i].v >= 0. ){ 
      phi_rad = acos(stars[i].u);
    }
    else{
      phi_rad = 2.*PI - acos(stars[i].u);
    }

    // calculate number of photons
    
    // mag_nphot() translates the star magnitude into number of photons,
    // using the expression log(flux)=-0.4*m-22.42 for each waveband
    // the resulting numbers ar stored in the array nph
    
    stars[i].mag_nphot(nph, pars.integtime_s, pars.mirr_radius_m, pars.verbose);
    
    // loop over all photons

    photinside=0; 

    for(k=0; k < 4; k++){ // loop over wavebands
    
      for(j=0; j<nph[k]; j++){ // loop over photons of this waveband

        // Check if we have overflowed the alllowed ph number

          if(totalphotinside >= iMAXPHOT){
            cout << "Warning: photon memory full. Can only store " << iMAXPHOT << " photons.\n";
            break;
            //exit(1);
          }
      
        // for every photon, a pair of uniform random x,y coordinates ( x*x+y*y<300 ) is generated
        // and a uniform random arrival time inside a time window given by the integration time.
            
        xdum_m=rand_coord(pars.mirr_radius_m);
        ydum_m=rand_coord(pars.mirr_radius_m);
      
        if((xdum_m*xdum_m+ydum_m*ydum_m)<pars.mirr_radius_m*pars.mirr_radius_m){
          
          randtime = rand_time(pars.integtime_s);         
          lambda_nm = rand_lambda(lmin_nm[k], lmax_nm[k]); 
          
          //fill the photon fields by using the member functions defined in photon.hxx
          
          photons[totalphotinside].starnum = starnumber;
          photons[totalphotinside].arrtime_sec = randtime;
          photons[totalphotinside].x_m = xdum_m;
          photons[totalphotinside].y_m = ydum_m; 
          photons[totalphotinside].u = stars[i].u;
          photons[totalphotinside].v = stars[i].v;
          photons[totalphotinside].lambda_nm = lambda_nm;

          if(pars.verbose > 2)
            cout << "PH " << starnumber << " " << randtime << " " << xdum_m << " "
                 << ydum_m << " " << stars[i].u << " " << stars[i].v << " " << lambda_nm 
                 << "\n";
          
          photinside++;

          totalphotinside++;

          totalnumphot++;        
          
        }
        else{
          
          totalnumphot++;
          continue;

        } // end if
        
      } // end photon loop   
    } // end waveband loop

    stars[i].numphot = photinside;
    //if (i>81) break;
    if (pars.verbose) cout<<"Star number= "<<i<< " (catalog number " << 
                        starnumber << ") Number of photons accepted = "<< photinside<<endl;
    
  } // end star loop

  cout << "Total number of photons accepted = " <<  totalphotinside << "\n";

  convertcorsika( 
                 ((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);

  return 0;
  
} // end main