Ignore:
Timestamp:
02/06/03 11:59:10 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc

    r1681 r1745  
    3838ClassImp(MObservatoryLocation);
    3939
    40 MObservatoryLocation::MObservatoryLocation(const char *name, const char *title)
     40void MObservatoryLocation::Init(const char *name, const char *title)
    4141{
    4242    fName  = name  ? name  : "MObservatoryLocation";
    4343    fTitle = title ? title : "Storage container for coordinates of an observatory";   
    44     //    TH1F fHorizon=new TH1F;   
    45     fgDegToRad=2*TMath::Pi()/360;   
    46     fLatitude  =  28.7594 * fgDegToRad; // rad; 28 45 34
    47     fLongitude =  17.8761 * fgDegToRad; // rad; 17 52 34;
     44}
     45
     46MObservatoryLocation::MObservatoryLocation(const char *name, const char *title)
     47{
     48    Init();
     49
     50    fLatitude  =  28.7594 / kRad2Deg; // rad; 28 45 34
     51    fLongitude =  17.8761 / kRad2Deg; // rad; 17 52 34;
    4852                                        // slalib uses + for WEST !!!
    4953    fElevation = 2300; // m
    5054    fObsName = "Observatorio del Roque de los Muchachos";
     55}
     56
     57MObservatoryLocation::MObservatoryLocation(LocationName_t name, const char *name=NULL, const char *title=NULL)
     58{
     59    Init();
     60
     61    switch (name)
     62    {
     63    case kMagic1:
     64    case kMagic2:
     65    case kRobertGarten:
     66        fLatitude  =  28.7594 / kRad2Deg; // rad; 28 45 34
     67        fLongitude =  17.8761 / kRad2Deg; // rad; 17 52 34;
     68                                          // slalib uses + for WEST !!!
     69        fElevation = 2300; // m
     70        fObsName = "Observatorio del Roque de los Muchachos";
     71        break;
     72    }
    5173}
    5274
    … …  
    6486  *fLog << all;
    6587  *fLog << fObsName << endl;
    66   *fLog << "Latitude " << (fLatitude > 0 ? (fLatitude/fgDegToRad) : -(fLatitude/fgDegToRad)) << " deg " << (fLatitude > 0 ? "W" : "E") << endl;
    67   *fLog << "Longitude " << (fLongitude > 0 ? (fLongitude/fgDegToRad) : -(fLongitude/fgDegToRad)) <<" deg " << (fLongitude < 0 ? "N" : "S") << endl;
     88  *fLog << "Latitude " << (fLatitude > 0 ? (fLatitude*kRad2Deg) : -(fLatitude*kRad2Deg)) << " deg " << (fLatitude > 0 ? "W" : "E") << endl;
     89  *fLog << "Longitude " << (fLongitude > 0 ? (fLongitude*kRad2Deg) : -(fLongitude*kRad2Deg)) <<" deg " << (fLongitude < 0 ? "N" : "S") << endl;
    6890  *fLog << "Elevation " << fElevation << "m" << endl;
    6991}
     92

  • trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h

    r1681 r1745  
    99{
    1010private:
    11   char* fObsName;
     11  TString fObsName;
    1212  Double_t fLatitude;
    1313  Double_t fLongitude;
    1414  Double_t fElevation;
    15   Double_t fgDegToRad;
     15  static Double_t fgDegToRad;
    1616  //  TH1F     fHorizon;
    1717
     18  void Init(const char *name, const char *title);
     19
    1820public:
     21    enum LocationName_t
     22    {
     23        kMagic1,
     24        kMagic2,
     25        kRobertGarten
     26    };
     27
    1928  MObservatoryLocation(const char *name=NULL, const char *title=NULL);
     29  MObservatoryLocation(LocationName_t name, const char *name=NULL, const char *title=NULL);
    2030  ~MObservatoryLocation();
    2131
    22   inline void SetLatitude(Double_t latitude) { fLatitude = latitude; }
    23   inline void SetLongitude(Double_t longitude) { fLongitude = longitude; }
    24   inline void SetElevation(Double_t elevation) { fElevation = elevation; }
    25   inline void SetObservatoryName(char* name) { fObsName = name; }
     32  void SetLatitude(Double_t latitude)  { fLatitude = latitude; }
     33  void SetLongitude(Double_t longitude) { fLongitude = longitude; }
     34  void SetElevation(Double_t elevation) { fElevation = elevation; }
     35  void SetObservatoryName(TString name) { fObsName = name; }
    2636 
    27   void MObservatoryLocation::Print(Option_t *) const;
     37  void Print(Option_t *) const;
    2838 
    29   inline Double_t GetLatitude() { return fLatitude/fgDegToRad; }
    30   inline Double_t GetLongitude() { return fLongitude/fgDegToRad; }
    31   inline Double_t GetElevation() { return fElevation; }
    32   Double_t GetLatitudeRad() { return fLatitude; }
    33   Double_t GetLongitudeRad() { return fLongitude; }
    34   char* GetObservatoryName() { return fObsName; }
     39  Double_t GetLatitude() const        { return fLatitude*kRad2Deg; }
     40  Double_t GetLongitude() const       { return fLongitude*kRad2Deg; }
     41  Double_t GetElevation() const      { return fElevation; }
     42  Double_t GetLatitudeRad() const     { return fLatitude; }
     43  Double_t GetLongitudeRad() const    { return fLongitude; }
     44  TString  GetObservatoryName() const { return fObsName; }
    3545  // Double_t GetHorizon(Double_t phi);
    3646  // void SetHorizonLine(TF1 hor) { fHorizon = hor; }

  • trunk/MagicSoft/Mars/mtemp/MVPObject.cc

    r1680 r1745  
    113113// HHMMsDDT, where RA is given in hours and minutes and Declination is
    114114// given by degrees DD and tenths of degrees T. "s" may be "+" or
    115 // "-"
    116 //
    117 void MVPObject::SetObjectByName(char* object)
     115// "-", eg. "1959+650"
     116//
     117void MVPObject::SetObjectByName(const char *object)
    118118{
    119119  fObjectName=object;
    120120  fGotName=kTRUE;
    121121 
    122 //   cout<<"OBJ:"<<object<<endl;
    123 
    124   unsigned int delim=0;
    125   for (unsigned int i=0; i<strlen(object); i++)
    126     if ((object[i]=='+')||(object[i]=='-'))
    127       delim=i;
    128  
    129   char ra[6];
    130   char de[6];
    131  
    132   unsigned int i;
    133   for (i=0;  i<=delim; i++)
    134     ra[i]=object[i];
    135   ra[i-1]=0;
    136 
    137   for (i=delim+1;  i<strlen(object); i++)
    138     de[i-delim-1]=object[i];
    139   de[i-delim-1]=0;
    140 
    141   Float_t RA, Dec;
    142 
    143   sscanf(ra,"%f",&RA);
    144   sscanf(de,"%f",&Dec);
    145 
    146 //   cout<<"OBJd:"<<Dec<<endl; //220
    147 //   cout<<"OBJr:"<<RA<<endl; //1959
    148 
    149   if (object[delim]=='-') Dec*=-1;
    150 
    151   fRA=(Double_t)( fgHrsToRad*  ((Int_t)(RA/100) + ( RA-(Int_t)(RA/100)*100)/60        ));
    152   fDec=(Double_t)( fgDegToRad* ((Int_t)(Dec/10) + (Dec-(Int_t)(Dec/10)*10 )/10        ));
    153 
    154  //  fRA=(Double_t)( fgHrsToRad*  ((Int_t)(RA/100)   + ((RA / 100)-(Int_t)(RA/100))/60   ));
    155 //   fDec=(Double_t)( fgDegToRad* ((Int_t)(Dec/10)  + ((Dec / 10)-(Int_t)(Dec/100))/10 ));
    156 
    157 //     cout<<"OBJd:"<<fDec/fgDegToRad<<endl;
    158 //     cout<<"OBJr:"<<fRA/fgHrsToRad<<endl;
     122  Int_t ra, dec;
     123  sscanf(object, "%d%d", &ra, &dec);
     124
     125  fRA  = fgHrsToRad * (0.01*ra + (r%100)/60.);
     126  fDec = fgDegToRad * 0.1 * dec;
    159127
    160128  fGotRA=kTRUE;

  • trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc

    r1681 r1745  
    190190  // each month [1]..[12]
    191191  Float_t visibility[13][18];
     192  memset(visibility, 0, 13*18*sizeof(Float_t));
     193  /*
    192194  for (int m=0;m<13;m++)
    193195    for (int z=0;z<18;z++)
    194196      visibility[m][z]=0;
    195 
     197  */
    196198  int fday, ftime;
    197199  Double_t phase=0;     
    … …  
    201203    Double_t todaysPhase=0;
    202204    Double_t moonIntensity;
    203     UInt_t i;
    204205    obsHours=0;
    205     for (i=0; i<daySlices; i++)
     206    for (UInt_t i=0; i<daySlices; i++)
    206207      { 
    207208        // Rearrange filling of bins such that a "day" doesn't start at midnight,
    … …  
    251252
    252253          // If moon is not up (or we should not use moon information...)
    253           if ((!fUseMoon)||(fMoon->GetAltitudeDeg()<=0)||(moonIntensity<60)) {
     254          if (!fUseMoon || fMoon->GetAltitudeDeg()<=0 || moonIntensity<60) {
    254255            // Fill MJD-UTC histogram
    255256            fMjdUtcYear->Fill(fgSecPerDay*(fday-fgMJD010170),fgSecPerDay*ftime/daySlices,fObject->GetAltitudeDeg());
    … …  
    296297
    297298      Double_t distance = fVenus->GetDistance(fObject);
    298       distance = distance < fMars->GetDistance(fObject) ? distance : fMars->GetDistance(fObject);
    299       distance = distance < fJupiter->GetDistance(fObject) ? distance : fJupiter->GetDistance(fObject);
    300       distance = distance < fSaturn->GetDistance(fObject) ? distance : fSaturn->GetDistance(fObject);
     299      distance = TMath::Min(distance, fMars->GetDistance(fObject));
     300      distance = TMath::Min(distance, fJupiter->GetDistance(fObject));
     301      distance = TMath::Min(distance, fSaturn->GetDistance(fObject));
    301302
    302303      fMjdPlanetDistance->Fill(fgSecPerDay*(fday-fgMJD010170),distance);     
    … …  
    319320      printf("%5d ",(Int_t)(visibility[m][z]/60));     
    320321    }
    321     cout<<endl;
     322    printf("\n");
    322323  }
    323324 
    324325  int vistimestart=0;
    325326  int vistimeend=0; 
    326   for (int m=1;m<13;m++) {
     327  for (int m=1; m<13; m++) {
    327328    int n = (m==1 ? 12 : m-1);
    328     if (((visibility[m][8]/60)>20)&&((visibility[n][8])/60 <= 20)) {
     329    if (visibility[m][8]/60>20 && visibility[n][8]/60<=20) {
    329330      vistimestart=m; // we're on the rising slope
    330331    }
    331332  }
    332333 
    333   for (int m=1;m<13;m++) {
     334  for (int m=1; m<13; m++) {
    334335    int n = (m==1 ? 12 : m-1);
    335     if (((visibility[m][8]/60)<20)&&((visibility[n][8]/60)>=20)) {
     336    if (visibility[m][8]/60<20 && visibility[n][8]/60>=20) {
    336337      vistimeend=n; // we're on the falling slope
    337338    }
    … …  
    341342
    342343
    343   gROOT->Reset(); 
     344  /*!!!!!!!!!!!!!!!!!!!!!!!*/gROOT->Reset();
    344345
    345346  TCanvas *cMjdUtcYear = new TCanvas("cMjdUtcYear", "Object Visibility MjdUtcYear", 1100,500);
    … …  
    430431  weekCounter++;
    431432 
    432   Int_t startday2 = fTime->MJDStartOfYear(year);
    433   Int_t stopday2  = fTime->MJDStartOfYear(year+1)-1;
    434  
    435     for (int day=startday2; day<stopday2+1; day=day+7)       
    436     {     
    437       simpleCounter[weekCounter]=weekCounter-1;
    438       objectHeight[weekCounter]=0;     
    439       for (int i=0; i<daySlices; i++)
    440         {       
    441           if (i>=48) {
    442             fday=day;
    443             ftime=i-48;
    444           } else {
    445             fday=day-1;
    446             ftime=i+48;
    447           }
    448          
    449           fTime->SetMJD(day,(Double_t)i/daySlices);       
    450           fObject->Process();                               
    451           fSun->Process();
    452          
    453           if (fSun->GetAltitudeDeg() < -25.0) {
    454             if (objectHeight[weekCounter]<(fObject->GetAltitudeDeg()))
    455               objectHeight[weekCounter]=fObject->GetAltitudeDeg();     
    456           }   
    457          
    458         } //i
    459       weekCounter++;
    460     } //day   
    461   simpleCounter[weekCounter]=weekCounter-2;
    462   objectHeight[weekCounter]=0;
    463   weekCounter++;
    464  
    465   TString months[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
    466  
    467   TCanvas *c2 = new TCanvas("c2", "Object Visibility",600,100);
     433  const Int_t startday2 = fTime->MJDStartOfYear(year);
     434  const Int_t stopday2  = fTime->MJDStartOfYear(year+1)-1;
     435 
     436    for (int day=startday2; day<stopday2+1; day+=7)
     437    {
     438        simpleCounter[weekCounter]=weekCounter-1;
     439        objectHeight[weekCounter]=0;
     440        for (int i=0; i<daySlices; i++)
     441        {
     442            if (i>=48)
     443            {
     444                fday=day;
     445                ftime=i-48;
     446            }
     447            else
     448            {
     449                fday=day-1;
     450                ftime=i+48;
     451            }
     452
     453            fTime->SetMJD(day,(Double_t)i/daySlices);
     454            fObject->Process();
     455            fSun->Process();
     456
     457            if (fSun->GetAltitudeDeg() < -25.0)
     458            {
     459                if (objectHeight[weekCounter]<(fObject->GetAltitudeDeg()))
     460                    objectHeight[weekCounter]=fObject->GetAltitudeDeg();
     461            }
     462
     463        } //i
     464        weekCounter++;
     465    } //day
     466    simpleCounter[weekCounter]=weekCounter-2;
     467    objectHeight[weekCounter]=0;
     468    weekCounter++;
     469
     470    TString months[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
     471
     472    TCanvas *c2 = new TCanvas("c2", "Object Visibility",600,100);
    468473 
    469   //  gStyle->SetOptTitle(0);
    470 //   gStyle->SetPadLeftMargin(0.000001);
    471 //   gStyle->SetPadRightMargin(0.000001);
    472 //   gStyle->SetPadTopMargin(0.001);
    473 //   gStyle->SetPadBottomMargin(0.000001);
    474   gPad->SetGrid();
    475  
    476   c2->SetGrid();
    477   TGraph *tg=new TGraph(weekCounter,simpleCounter,objectHeight);
    478  
    479   tg->SetMinimum(1);
    480   tg->SetMaximum(90);
    481  
    482   Double_t maxza = abs(fObject->GetDec()-fObs->GetLatitude());
    483  
    484   if (maxza > 90) maxza=90;
    485   if (maxza < 0) maxza=00;
    486  
    487   cout << "MaxZA:    ";
    488   cout << maxza << endl;
    489  
    490   if (maxza < 30) { //colors=green to yellow
    491     maxza=9*maxza/30;
    492     maxza=80+maxza;
    493    
    494   } else { //colors=yellow to red
    495     maxza-=30;
    496     maxza=11*maxza/60;
    497     maxza=89+maxza;
    498   }   
    499  
     474    //  gStyle->SetOptTitle(0);
     475    //   gStyle->SetPadLeftMargin(0.000001);
     476    //   gStyle->SetPadRightMargin(0.000001);
     477    //   gStyle->SetPadTopMargin(0.001);
     478    //   gStyle->SetPadBottomMargin(0.000001);
     479    gPad->SetGrid();
     480
     481    c2->SetGrid();
     482    TGraph *tg=new TGraph(weekCounter,simpleCounter,objectHeight);
     483
     484    tg->SetMinimum(1);
     485    tg->SetMaximum(90);
     486
     487    Double_t maxza = abs(fObject->GetDec()-fObs->GetLatitude());
     488
     489    if (maxza > 90) maxza=90;
     490    if (maxza < 0) maxza=00;
     491
     492    cout << "MaxZA:    ";
     493    cout << maxza << endl;
     494
     495    if (maxza < 30) { //colors=green to yellow
     496        maxza *= 9/30;
     497        maxza += 80;
     498
     499    } else { //colors=yellow to red
     500        maxza -= 30;
     501        maxza *= 11/60;
     502        maxza += 89;
     503    }
     504
    500505  tg->SetFillColor((Int_t)maxza);
    501506  tg->SetLineColor((Int_t)maxza);
    … …  
    515520 
    516521  if ((vistimestart<13)&&(vistimeend>0)) {
    517     TString label;
    518     label=months[vistimestart-1]+"-"+months[vistimeend-1];
     522    TString label=months[vistimestart-1]+"-"+months[vistimeend-1];
    519523    TPaveLabel* l2= new TPaveLabel(35,46,50,88, label);
    520524    l2->SetBorderSize(0);
    … …  
    528532  c2->Update();
    529533 
    530 return kTRUE;
     534  return kTRUE;
    531535}
    532536
    … …  
    553557    istar =  pow(10.,istar);
    554558    if(fabs(alpha) < 7.)   /* crude accounting for opposition effect */
    555         istar = istar * (1.35 - 0.05 * fabs(istar));
     559        istar *= 1.35 - 0.05 * fabs(istar);
    556560        /* 35 per cent brighter at full, effect tapering linearly to
    557561           zero at 7 degrees away from full. mentioned peripherally in
    … …  
    559563    double fofrho = 229087. * (1.06 + cos(rho_rad)*cos(rho_rad));
    560564    if(fabs(rho) > 10.)
    561        fofrho=fofrho+pow(10.,(6.15 - rho/40.));            /* eqn 21 */
     565        fofrho+=pow(10.,(6.15 - rho/40.));            /* eqn 21 */
    562566    else if (fabs(rho) > 0.25)
    563        fofrho= fofrho+ 6.2e7 / (rho*rho);   /* eqn 19 */
     567        fofrho+=6.2e7 / (rho*rho);   /* eqn 19 */
    564568    else fofrho = fofrho+9.9e8;  /*for 1/4 degree -- radius of moon! */
     569
    565570    double Xzm = sqrt(1.0 - 0.96*sin(Zmoon)*sin(Zmoon));
     571
    566572    if(Xzm != 0.) Xzm = 1./Xzm;
    567           else Xzm = 10000.;
     573    else Xzm = 10000.;
     574
    568575    double Xo = sqrt(1.0 - 0.96*sin(Z)*sin(Z));
    569576    if(Xo != 0.) Xo = 1./Xo;
    570           else Xo = 10000.;
     577    else Xo = 10000.;
     578
    571579    double Bmoon = fofrho * istar * pow(10.,(-0.4*kzen*Xzm))
    572           * (1. - pow(10.,(-0.4*kzen*Xo)));   /* nanoLamberts */
     580        * (1. - pow(10.,(-0.4*kzen*Xo)));   /* nanoLamberts */
    573581    //    cout << " Bmoon=" << Bmoon;
    574582    if(Bmoon > 0.001)

  • trunk/MagicSoft/Mars/mtemp/MVPPlotter.h

    r1681 r1745  
    11#ifndef MARS_MVPPlotter
    22#define MARS_MVPPlotter
    3 
    4 #ifndef MARS_MTask
    5 #include "MTask.h"
    6 #endif
    7 
    8 #ifndef MARS_MParList
    9 #include "MParList.h"
    10 #endif
    113
    124#ifndef MARS_MTaskList
Note: See TracChangeset for help on using the changeset viewer.