/[PAMELA software]/DarthVader/OrbitalInfo/src/OrbitalInfoCore.cpp
ViewVC logotype

Diff of /DarthVader/OrbitalInfo/src/OrbitalInfoCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.77 by mocchiutti, Tue Sep 16 19:47:34 2014 UTC revision 1.84 by mocchiut, Thu Mar 26 14:55:38 2015 UTC
# Line 26  Line 26 
26  #include <TSQLServer.h>  #include <TSQLServer.h>
27  #include <TSQLRow.h>  #include <TSQLRow.h>
28  #include <TSQLResult.h>  #include <TSQLResult.h>
29    #include <TObjectTable.h>
30  //  //
31  // RunInfo header  // RunInfo header
32  //  //
# Line 61  using namespace std; Line 62  using namespace std;
62  // CORE ROUTINE  // CORE ROUTINE
63  //  //
64  //  //
65  int OrbitalInfoCore(UInt_t run, TFile *file, GL_TABLES *glt, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){  int OrbitalInfoCore(UInt_t run, TFile *file, GL_TABLES *glt, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){  
66    //    //
67    Int_t i = 0;    Int_t i = 0;
68    TString host = glt->CGetHost();    TString host = glt->CGetHost();
# Line 71  int OrbitalInfoCore(UInt_t run, TFile *f Line 72  int OrbitalInfoCore(UInt_t run, TFile *f
72    //    //
73    stringstream myquery;    stringstream myquery;
74    myquery.str("");    myquery.str("");
75    myquery << "SET time_zone='+0:00'";    myquery << "SET time_zone='+0:00';";
76    delete dbc->Query(myquery.str().c_str());    delete dbc->Query(myquery.str().c_str());
77      delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
78    //    //
79    TString processFolder = Form("OrbitalInfoFolder_%u",run);    TString processFolder = Form("OrbitalInfoFolder_%u",run);
80    //    //
# Line 90  int OrbitalInfoCore(UInt_t run, TFile *f Line 92  int OrbitalInfoCore(UInt_t run, TFile *f
92        if ( !strcmp(OrbitalInfoargv[i],"-processFolder") ) {        if ( !strcmp(OrbitalInfoargv[i],"-processFolder") ) {
93          if ( OrbitalInfoargc < i+1 ){          if ( OrbitalInfoargc < i+1 ){
94            throw -3;            throw -3;
95          };          }
96          processFolder = (TString)OrbitalInfoargv[i+1];          processFolder = (TString)OrbitalInfoargv[i+1];
97          i++;          i++;
98        };        }
99        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {
100          verbose = true;          verbose = true;
101          debug = true;          debug = true;
102        };        }
103        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
104          verbose = true;          verbose = true;
105        };        }
106        if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {
107          standalone = true;          standalone = true;
108        };        }
109        if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {
110          standalone = false;          standalone = false;
111        };        }
112        i++;        i++;
113      };      }
114    };    }
115      if ( debug ){
116        printf("START\n");
117        gObjectTable->Print();
118      }
119    //    //
120    const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));    const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));
121    //    //
# Line 189  int OrbitalInfoCore(UInt_t run, TFile *f Line 195  int OrbitalInfoCore(UInt_t run, TFile *f
195    //    //
196    // IGRF stuff    // IGRF stuff
197    //    //
198    Double_t dimo = 0.0; // dipole moment (computed from dat files) // EM GCC 4.7    Float_t dimo = 0.0; // dipole moment (computed from dat files) // EM GCC 4.7
199    Float_t bnorth, beast, bdown, babs;    Float_t bnorth, beast, bdown, babs;
200    Float_t xl; // L value    Float_t xl; // L value
201    Float_t icode; // code value for L accuracy (see fortran code)    Float_t icode; // code value for L accuracy (see fortran code)
# Line 242  int OrbitalInfoCore(UInt_t run, TFile *f Line 248  int OrbitalInfoCore(UInt_t run, TFile *f
248    // Initialize fortran routines!!!    // Initialize fortran routines!!!
249    Int_t ltp1 = 0;    Int_t ltp1 = 0;
250    Int_t ltp2 = 0;    Int_t ltp2 = 0;
   Int_t ltp3 = 0;  
   //  Int_t uno = 1;  
   //  const char *niente = " ";  
251    GL_PARAM *glparam0 = new GL_PARAM();    GL_PARAM *glparam0 = new GL_PARAM();
252    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
253    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
   GL_PARAM *glparam3 = new GL_PARAM();  
254    
255    //    //
256    // Orientation variables. Vitaly    // Orientation variables. Vitaly
# Line 281  int OrbitalInfoCore(UInt_t run, TFile *f Line 283  int OrbitalInfoCore(UInt_t run, TFile *f
283    vector<Float_t> recq2;    vector<Float_t> recq2;
284    vector<Float_t> recq3;    vector<Float_t> recq3;
285    Float_t Norm = 1;    Float_t Norm = 1;
286      recqtime.reserve(1500000);
287      recq0.reserve(1500000);
288      recq1.reserve(1500000);
289      recq2.reserve(1500000);
290      recq3.reserve(1500000);
291    
292    vector<UInt_t> RTtime1;    vector<UInt_t> RTtime1;
293    vector<UInt_t> RTtime2;    vector<UInt_t> RTtime2;
# Line 293  int OrbitalInfoCore(UInt_t run, TFile *f Line 300  int OrbitalInfoCore(UInt_t run, TFile *f
300    vector<UInt_t> RTpluto2;    vector<UInt_t> RTpluto2;
301    vector<UInt_t> RTpluto1;    vector<UInt_t> RTpluto1;
302    vector<Int_t> RTerrq;    vector<Int_t> RTerrq;
303      vector<Int_t> RTqual;
304      RTtime1.reserve(200000);
305      RTtime2.reserve(200000);
306      RTbank1.reserve(200000);
307      RTbank2.reserve(200000);
308      RTbpluto1.reserve(200000);
309      RTbpluto2.reserve(200000);
310      RTazim.reserve(200000);
311      RTstart.reserve(200000);
312      RTpluto1.reserve(200000);
313      RTpluto2.reserve(200000);
314      RTerrq.reserve(200000);
315      RTqual.reserve(200000);
316    
317    TClonesArray *tcNucleiTrk = NULL;    TClonesArray *tcNucleiTrk = NULL;
318    TClonesArray *tcExtNucleiTrk = NULL;    TClonesArray *tcExtNucleiTrk = NULL;
# Line 336  int OrbitalInfoCore(UInt_t run, TFile *f Line 356  int OrbitalInfoCore(UInt_t run, TFile *f
356      in>>recq2[sizee-1];      in>>recq2[sizee-1];
357      in>>recq3[sizee-1];      in>>recq3[sizee-1];
358      in>>Norm;      in>>Norm;
359    /* CHECK RECOVERED QUATERNIONS PROBLEM
360        if(recqtime[sizee-1]>=1160987921.75 && recqtime[sizee-1]<=1160987932.00){
361          cout<<"We found it at start"<<"\t"<<recqtime[sizee-1]<<endl;
362        } */
363    }    }
364    in.close();    in.close();
365    if ( verbose ) cout<<"We have read recovered data"<<endl;    if ( verbose ) cout<<"We have read recovered data"<<endl;
366    if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;    if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;
367      if ( debug )  printf(" RQ size %i RQ capacity %i \n",(int)recqtime.size(),(int)recqtime.capacity());
368      
369    if ( verbose ) cout<<"read Rotation Table"<<endl;    if ( verbose ) cout<<"read Rotation Table"<<endl;
370        
371    parerror2=glparam0->Query_GL_PARAM(1,305,dbc);    parerror2=glparam0->Query_GL_PARAM(1,305,dbc);
# Line 360  int OrbitalInfoCore(UInt_t run, TFile *f Line 385  int OrbitalInfoCore(UInt_t run, TFile *f
385      RTstart.resize(sizee+1);      RTstart.resize(sizee+1);
386      RTpluto1.resize(sizee+1);      RTpluto1.resize(sizee+1);
387      RTbpluto1.resize(sizee+1);      RTbpluto1.resize(sizee+1);
388        RTqual.resize(sizee+1);
389      an>>RTtime1[sizee-1];      an>>RTtime1[sizee-1];
390      an>>RTbank1[sizee-1];      an>>RTbank1[sizee-1];
391      an>>RTstart[sizee-1];      an>>RTstart[sizee-1];
# Line 367  int OrbitalInfoCore(UInt_t run, TFile *f Line 393  int OrbitalInfoCore(UInt_t run, TFile *f
393      an>>RTbpluto1[sizee-1];      an>>RTbpluto1[sizee-1];
394      an>>RTazim[sizee-1];      an>>RTazim[sizee-1];
395      an>>RTerrq[sizee-1];      an>>RTerrq[sizee-1];
396        an>>RTqual[sizee-1];
397      if(sizee>1) {      if(sizee>1) {
398        RTtime2.resize(sizee+1);        RTtime2.resize(sizee+1);
399        RTbank2.resize(sizee+1);        RTbank2.resize(sizee+1);
# Line 385  int OrbitalInfoCore(UInt_t run, TFile *f Line 412  int OrbitalInfoCore(UInt_t run, TFile *f
412        
413    if ( verbose ) cout<<"We have read Rotation Table"<<endl;    if ( verbose ) cout<<"We have read Rotation Table"<<endl;
414      //Geomagnetic coordinates calculations staff      //Geomagnetic coordinates calculations staff
415      
416      if ( debug ) printf(" RT size %i RT capacity %i \n",(int)RTtime2.size(),(int)RTtime2.capacity());
417    
418    GMtype_CoordGeodetic location;    GMtype_CoordGeodetic location;
419    //  GMtype_CoordDipole GMlocation;    //  GMtype_CoordDipole GMlocation;
# Line 395  int OrbitalInfoCore(UInt_t run, TFile *f Line 424  int OrbitalInfoCore(UInt_t run, TFile *f
424    //    TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";    //    TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
425    //  }    //  }
426    
   //  GM_ScanIGRF(glparam->PATH, &G0, &G1, &H1);  
   GM_ScanIGRF(dbc, &G0, &G1, &H1);  
   
427    //cout << G0.element[0] << "\t" << G1.element[0] << "\t" << H1.element[0] << endl;    //cout << G0.element[0] << "\t" << G1.element[0] << "\t" << H1.element[0] << endl;
428    //cout << G0.element[5] << "\t" << G1.element[5] << "\t" << H1.element[5] << endl;    //cout << G0.element[5] << "\t" << G1.element[5] << "\t" << H1.element[5] << endl;
429    
# Line 707  int OrbitalInfoCore(UInt_t run, TFile *f Line 733  int OrbitalInfoCore(UInt_t run, TFile *f
733    // Get the list of run to be processed, if only one run has to be processed the list will contain one entry only.    // Get the list of run to be processed, if only one run has to be processed the list will contain one entry only.
734    //    //
735    runlist = runinfo->GetRunList();    runlist = runinfo->GetRunList();
736      if ( debug ){
737        printf("BEFORE LOOP ON RUN\n");
738        gObjectTable->Print();
739      }
740    //    //
741    // Loop over the run to be processed    // Loop over the run to be processed
742    //    //
743    for (UInt_t irun=0; irun < numbofrun; irun++){    for (UInt_t irun=0; irun < numbofrun; irun++){ //===>
744    
745      L_QQ_Q_l_lower = new Quaternions();      L_QQ_Q_l_lower = new Quaternions();
746      RYPang_lower = new InclinationInfo();      RYPang_lower = new InclinationInfo();
# Line 819  int OrbitalInfoCore(UInt_t run, TFile *f Line 849  int OrbitalInfoCore(UInt_t run, TFile *f
849        goto closeandexit;        goto closeandexit;
850      };      };
851    
     //  
     // open IGRF files and do it only once if we are processing a full level2 file  
     //  
     if ( !igrfloaded ){  
   
       if ( l0head->GetEntry(runinfo->EV_FROM) > 0 ){  
         igrfloaded = true;  
         //  
         // absolute time of first event of the run (it should not matter a lot)  
         //  
         ph = eh->GetPscuHeader();  
         atime = dbtime->DBabsTime(ph->GetOrbitalTime());  
           
         parerror=glparam->Query_GL_PARAM(atime-anni5,301,dbc); // parameters stored in DB in GL_PRAM table    
         if ( parerror<0 ) {  
           code = parerror;  
           goto closeandexit;  
         }  
         ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();  
         if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());  
         //  
         parerror=glparam2->Query_GL_PARAM(atime,301,dbc); // parameters stored in DB in GL_PRAM table    
         if ( parerror<0 ) {  
           code = parerror;  
           goto closeandexit;  
         }  
         ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();  
         if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());  
         //  
         parerror=glparam3->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table  
         if ( parerror<0 ) {  
           code = parerror;  
           goto closeandexit;  
         }  
         ltp3 = (Int_t)(glparam3->PATH+glparam3->NAME).Length();  
         if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam3->PATH+glparam3->NAME).Data());  
         //  
         initize_((char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2,(char *)(glparam3->PATH+glparam3->NAME).Data(),&ltp3);  
         //  
         if (debug) cout<<"initize: "<<(char *)(glparam->PATH+glparam->NAME).Data()<<"\t"<<(char *)(glparam2->PATH+glparam2->NAME).Data()<<"\t"<<(char *)(glparam3->PATH+glparam3->NAME).Data()<<endl;  
       }  
     }  
     //  
     // End IGRF stuff//  
     //  
   
852      ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();      ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
853      ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);      ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
854      ULong_t DeltaOBT = TimeSync - ObtSync;      ULong_t DeltaOBT = TimeSync - ObtSync;
# Line 900  int OrbitalInfoCore(UInt_t run, TFile *f Line 884  int OrbitalInfoCore(UInt_t run, TFile *f
884          //          //
885          l0fid[i] = (UInt_t)atoll(Row->GetField(0));          l0fid[i] = (UInt_t)atoll(Row->GetField(0));
886          i--;          i--;
887            if (Row){  // memleak!
888              delete Row;
889              Row = 0;
890            }
891          Row = pResult->Next();            Row = pResult->Next();  
892          //          //
893        };        }
894          if (Row) delete Row;
895        pResult->Delete();        pResult->Delete();
896      };      }
897      //      //
898      myquery.str("");      myquery.str("");
899      myquery << "select ID_ROOT_L0 from GL_RUN where RUNHEADER_TIME>" << runinfo->RUNHEADER_TIME << " group by ID_ROOT_L0 order by RUNHEADER_TIME asc limit 5;";      myquery << "select ID_ROOT_L0 from GL_RUN where RUNHEADER_TIME>" << runinfo->RUNHEADER_TIME << " group by ID_ROOT_L0 order by RUNHEADER_TIME asc limit 5;";
# Line 922  int OrbitalInfoCore(UInt_t run, TFile *f Line 911  int OrbitalInfoCore(UInt_t run, TFile *f
911          //          //
912          l0fid[i] = (UInt_t)atoll(Row->GetField(0));          l0fid[i] = (UInt_t)atoll(Row->GetField(0));
913          i++;          i++;
914            if (Row){  // memleak!
915              delete Row;
916              Row = 0;
917            }
918          Row = pResult->Next();            Row = pResult->Next();  
919          //          //
920        };        }
921          if (Row) delete Row;
922        pResult->Delete();        pResult->Delete();
923      };      }
924      //      //
925      i = 0;      i = 0;
926      UInt_t previd = 0;      UInt_t previd = 0;
# Line 945  int OrbitalInfoCore(UInt_t run, TFile *f Line 939  int OrbitalInfoCore(UInt_t run, TFile *f
939            if ( debug ) printf(" Using inclination informations from file: %s \n",(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1)).Data());            if ( debug ) printf(" Using inclination informations from file: %s \n",(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1)).Data());
940            ch->Add(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1));            ch->Add(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1));
941            //            //
942              if (Row) delete Row;
943            pResult->Delete();            pResult->Delete();
944          };          }
945        };        }
946        i++;        i++;
947      };      }
948      //      //
949      ch->SetBranchAddress("Mcmd",&mcmdev);      ch->SetBranchAddress("Mcmd",&mcmdev);
950      neventsm = ch->GetEntries();      neventsm = ch->GetEntries();
# Line 976  int OrbitalInfoCore(UInt_t run, TFile *f Line 971  int OrbitalInfoCore(UInt_t run, TFile *f
971      vector<Float_t> qYaw;      vector<Float_t> qYaw;
972      vector<Int_t> qmode;      vector<Int_t> qmode;
973    
974        q0.reserve(4096);
975        q1.reserve(4096);
976        q2.reserve(4096);
977        q3.reserve(4096);
978        qtime.reserve(4096);
979        qPitch.reserve(4096);
980        qRoll.reserve(4096);
981        qYaw.reserve(4096);
982        qmode.reserve(4096);
983        if ( debug ) printf(" q0 capa %i \n",(int)q0.capacity());
984      Int_t nt = 0;      Int_t nt = 0;
985      UInt_t must = 0;      UInt_t must = 0;
986    
987        Int_t currentYear = 0;
988        Int_t previousYear = 0;
989    
990      //      //
991      // run over all the events of the run      // run over all the events of the run
992      //      //
993      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");
994        if ( debug ){
995          printf("BEFORE LOOP ON EVENTS\n");
996          gObjectTable->Print();
997        }
998      //      //
999      //      //
1000      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
1001          //for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+10); re++){
1002    
1003        //        //
1004        if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);          if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
1005        if ( debug ) printf(" %i \n",procev);              if ( debug ) printf(" %i \n",procev);      
# Line 1007  int OrbitalInfoCore(UInt_t run, TFile *f Line 1021  int OrbitalInfoCore(UInt_t run, TFile *f
1021          continue;          continue;
1022        }        }
1023    
1024          // just for testing:
1025          //      if (re >= 5+runinfo->EV_FROM) atime += anni5;
1026          //      if (re >= 7+runinfo->EV_FROM) atime += anni5;
1027          //      if (re >= 9+runinfo->EV_FROM) atime += anni5;
1028    
1029          //
1030          // open IGRF files and do it only once if we are processing a full level2 file
1031          //
1032          Float_t kkyear;
1033          UInt_t kyear, kmonth, kday, khour, kmin, ksec;
1034          //
1035          TTimeStamp kt = TTimeStamp(atime, kTRUE);
1036          kt.GetDate(kTRUE, 0, &kyear, &kmonth, &kday);
1037          kt.GetTime(kTRUE, 0, &khour, &kmin, &ksec);
1038          kkyear = (float) kyear
1039            + (kmonth*31.+ (float) kday)/365.
1040            + (khour*3600.+kmin*60.+(float)ksec)/(24.*3600.*365.);
1041          currentYear = int(kkyear/5.) * 5;
1042          if ( debug ) printf(" prevy %i curry %i igrfloaded %i \n",previousYear,currentYear,igrfloaded);
1043          if ( currentYear != previousYear ) igrfloaded = false;
1044          previousYear = currentYear;
1045          if ( debug ) printf(" prevy %i curry %i igrfloaded %i \n",previousYear,currentYear,igrfloaded);
1046          //
1047          if ( !igrfloaded ){
1048            
1049            igrfloaded = true;
1050            
1051            parerror=glparam->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table  
1052            if ( parerror<0 ) {
1053              code = parerror;
1054              goto closeandexit;
1055            }
1056            ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();
1057            if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
1058            //
1059            if ( glparam->NAME.EndsWith("s.txt") || glparam->NAME.EndsWith("s.dat") ){
1060              if ( verbose ) printf("ERROR: Current date is beyond the latest secular variation file time span. Please update IGRF files to process data\n");
1061              throw -906;
1062            }
1063            //
1064            int isSecular = false;
1065            //
1066            parerror=glparam2->Query_GL_PARAM(atime+anni5,302,dbc); // parameters stored in DB in GL_PRAM table  
1067            if ( parerror<0 ) {
1068              code = parerror;
1069              goto closeandexit;
1070            }
1071            ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();
1072            if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());
1073            if ( glparam2->NAME.EndsWith("s.txt") || glparam2->NAME.EndsWith("s.dat") ){
1074              isSecular = true;
1075              if ( verbose ) printf(" Using secular variation file and hence fortran subroutine 'extrapolation'\n");
1076            } else {
1077              if ( verbose ) printf(" Using two field measurement files and hence fortran subroutine 'interpolation'\n");
1078            }
1079            //
1080            initize_(&isSecular,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2);
1081            //
1082            if (debug) cout<<"initize: "<<(char *)(glparam->PATH+glparam->NAME).Data()<<"\t"<<(char *)(glparam2->PATH+glparam2->NAME).Data()<<"\t isSecular? "<<isSecular<<endl;  
1083    
1084            //        GM_ScanIGRF(dbc, &G0, &G1, &H1);
1085            TString igrfFile1 = glparam->PATH+glparam->NAME;
1086            TString igrfFile2 = glparam2->PATH+glparam2->NAME;
1087            GM_SetIGRF(isSecular,igrfFile1,igrfFile2, &G0, &G1, &H1);
1088          }
1089          //
1090          // End IGRF stuff//
1091          //
1092    
1093        //        //
1094        // retrieve tof informations        // retrieve tof informations
1095        //        //
# Line 1122  int OrbitalInfoCore(UInt_t run, TFile *f Line 1205  int OrbitalInfoCore(UInt_t run, TFile *f
1205        cCoordGeo coo;        cCoordGeo coo;
1206        Float_t jyear=0.;            Float_t jyear=0.;    
1207        //        //
1208        if(atime >= gltle->GetToTime()) {        if(atime >= gltle->GetToTime() || atime < gltle->GetFromTime() ) {  // AGH! bug when reprocessing??
1209          if ( !gltle->Query(atime, dbc) ){          if ( !gltle->Query(atime, dbc) ){
1210            //                  //      
1211            // Compute the magnetic dipole moment.            // Compute the magnetic dipole moment.
# Line 1142  int OrbitalInfoCore(UInt_t run, TFile *f Line 1225  int OrbitalInfoCore(UInt_t run, TFile *f
1225            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
1226            if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);            if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);
1227    
1228            GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);            //      GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);
1229              GM_TimeAdjustCoefs(GM_STARTYEAR, (jyear-currentYear+GM_STARTYEAR), G0, G1, H1, &Model);  // EM: input this way due to the new way of storing data into Gn,H1 and to avoid changing GM_Time...
1230            GM_PoleLocation(Model, &Pole);            GM_PoleLocation(Model, &Pole);
1231                        
1232          } else {          } else {
# Line 1265  int OrbitalInfoCore(UInt_t run, TFile *f Line 1349  int OrbitalInfoCore(UInt_t run, TFile *f
1349                      for(Int_t mu = nt;mu<recSize;mu++){                      for(Int_t mu = nt;mu<recSize;mu++){
1350                        if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){                        if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1351                           if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){                           if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1352                             nt=mu;  //                         nt=mu;
1353                             Int_t sizeqmcmd = qtime.size();                             Int_t sizeqmcmd = qtime.size();
1354                             inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                             inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1355                             qtime[sizeqmcmd]=recqtime[mu];                             qtime[sizeqmcmd]=recqtime[mu];
# Line 1279  int OrbitalInfoCore(UInt_t run, TFile *f Line 1363  int OrbitalInfoCore(UInt_t run, TFile *f
1363                             qRoll[sizeqmcmd]=RYPang_upper->Kren;                             qRoll[sizeqmcmd]=RYPang_upper->Kren;
1364                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1365                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1366    /* CHECK RECOVERED QUATERNIONS PROBLEM */
1367    if(recqtime[mu]>=1160987921.75 && recqtime[mu]<=1160987932.00){
1368      cout<<"We found it while checking all quaternions"<<"\t"<<recqtime[mu]<<endl;
1369    }
1370                           }                           }
1371                        }                        }
1372                        if(recqtime[mu]>=u_time){                        if(recqtime[mu]>=u_time){
# Line 1480  int OrbitalInfoCore(UInt_t run, TFile *f Line 1568  int OrbitalInfoCore(UInt_t run, TFile *f
1568        // 10REDNEW        // 10REDNEW
1569        Int_t errq=0;        Int_t errq=0;
1570        Int_t azim=0;        Int_t azim=0;
1571          Int_t qual=0;
1572        Int_t MU=0;        Int_t MU=0;
1573        for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){        for(UInt_t mu = 0;mu<RTtime2.size()-1;mu++){
1574          if(atime<=RTstart[mu+1] && atime>=RTstart[mu]){          if(atime<RTstart[mu+1] && atime>=RTstart[mu]){
1575            errq=RTerrq[mu];            errq=RTerrq[mu];
1576            azim=RTazim[mu];            azim=RTazim[mu];
1577              qual=RTqual[mu];
1578            MU=mu;            MU=mu;
1579            break;            break;
1580          }          }
1581        }        }
1582        orbitalinfo->errq = errq;        orbitalinfo->errq = errq;
1583        orbitalinfo->azim = azim;        orbitalinfo->azim = azim;
1584          orbitalinfo->rtqual=qual;
1585        orbitalinfo->qkind = 0;        orbitalinfo->qkind = 0;
1586                
1587        if ( debug ) printf(" coord done \n");        if ( debug ) printf(" coord done \n");
# Line 1537  int OrbitalInfoCore(UInt_t run, TFile *f Line 1628  int OrbitalInfoCore(UInt_t run, TFile *f
1628          orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //          orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //
1629          if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff:  "<< orbitalinfo->cutoffsvl << endl;          if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff:  "<< orbitalinfo->cutoffsvl << endl;
1630    
1631          /*          
1632            ---------- Forwarded message ----------  //           ---------- Forwarded message ----------
1633            Date: Wed, 09 May 2012 12:16:47 +0200  //           Date: Wed, 09 May 2012 12:16:47 +0200
1634            From: Alessandro Bruno <alessandro.bruno@ba.infn.it>  //           From: Alessandro Bruno <alessandro.bruno@ba.infn.it>
1635            To: Mirko Boezio <mirko.boezio@ts.infn.it>  //           To: Mirko Boezio <mirko.boezio@ts.infn.it>
1636            Cc: Francesco S. Cafagna <Francesco.Cafagna@ba.infn.it>  //           Cc: Francesco S. Cafagna <Francesco.Cafagna@ba.infn.it>
1637            Subject: Störmer vertical cutoff  //           Subject: Störmer vertical cutoff
1638    
1639            Ciao Mirko,  //           Ciao Mirko,
1640            volevo segnalarti che il valore dello Störmer vertical cutoff nel Level2 è  //           volevo segnalarti che il valore dello Störmer vertical cutoff nel Level2 è
1641            sovrastimato di circa il 4%.  //           sovrastimato di circa il 4%.
1642            Dopo un'approfondita analisi con l'IGRF-05 abbiamo ricavano un valore pari  //           Dopo un'approfondita analisi con l'IGRF-05 abbiamo ricavano un valore pari
1643            a: 14.295 / L^2 anzichè 14.9 / L^2, valore obsoleto in quanto riferito agli  //           a: 14.295 / L^2 anzichè 14.9 / L^2, valore obsoleto in quanto riferito agli
1644            anni '50.  //           anni '50.
1645          */  //        
1646          //14.9/(xl*xl);          //14.9/(xl*xl);
1647          orbitalinfo->igrf_icode = icode;          orbitalinfo->igrf_icode = icode;
1648          //          //
# Line 1568  int OrbitalInfoCore(UInt_t run, TFile *f Line 1659  int OrbitalInfoCore(UInt_t run, TFile *f
1659          Float_t By = orbitalinfo->Beast;          Float_t By = orbitalinfo->Beast;
1660          Float_t Bz = orbitalinfo->Bnorth;          Float_t Bz = orbitalinfo->Bnorth;
1661    
1662          TMatrixD Qiji(3,3);  //      TMatrixD Qiji(3,3);
1663          TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);          TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);
1664          TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);          TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1665    
1666  //10REDNEW  //10REDNEW
1667          /* If initial orientation data have reason to be inaccurate */          // If initial orientation data have reason to be inaccurate
1668          Float_t tg = 0.00;         Float_t tg = 0.00;
1669         Float_t tmptg;         Float_t tmptg;
1670         if(MU!=0){         if(MU!=0){
1671  //      if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){               // 10RED CHECK  (comparison between three metod of recovering orientation)  //      if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){               // 10RED CHECK  (comparison between three metod of recovering orientation)
1672         if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && TMath::Abs(orbitalinfo->etha)>0.01) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || orbitalinfo->TimeGap>10.0 || ((modf(orbitalinfo->TimeGap,&tmptg)*1000>10 || modf(orbitalinfo->TimeGap,&tmptg)*1000==0.0) && orbitalinfo->TimeGap>2.0)))){         if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && TMath::Abs(orbitalinfo->etha)>0.1) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || orbitalinfo->TimeGap>10.0 || ((modf(orbitalinfo->TimeGap,&tmptg)*1000>10 || modf(orbitalinfo->TimeGap,&tmptg)*1000==0.0) && orbitalinfo->TimeGap>2.0)))){
1673          /*  found in Rotation Table this data for this time interval*/          //found in Rotation Table this data for this time interval
1674          if(atime<RTtime1[0])          if(atime<RTtime1[0])
1675            orbitalinfo->azim = 5;        //means that RotationTable no started yet            orbitalinfo->azim = 5;        //means that RotationTable no started yet
1676         else{         else{
# Line 1595  int OrbitalInfoCore(UInt_t run, TFile *f Line 1686  int OrbitalInfoCore(UInt_t run, TFile *f
1686                  bank=kar*atime+bak;                  bank=kar*atime+bak;
1687                }                }
1688                if(atime>=RTstart[MU] && atime<RTpluto1[MU]){                if(atime>=RTstart[MU] && atime<RTpluto1[MU]){
1689                   Double_t s_dBdt2=(RTbpluto1[MU]-RTbank1[MU])/(RTpluto1[MU]-RTstart[MU]);                   Double_t s_dBdt2=(RTbpluto1[MU]-RTbank1[MU])/(Int_t)(RTpluto1[MU]-RTstart[MU]);
1690                   Double_t s_t2=((Double_t)RTpluto1[MU]+(Double_t)RTstart[MU])/2. - RTstart[MU];                   Double_t s_t2=((Double_t)RTpluto1[MU]+(Double_t)RTstart[MU])/2. - RTstart[MU];
1691                   Double_t s_t1=RTstart[MU]-RTstart[MU];                   Double_t s_t1=RTstart[MU]-RTstart[MU];
1692                   Double_t s_k=s_dBdt2/(s_t2-s_t1);                   Double_t s_k=s_dBdt2/(s_t2-s_t1);
# Line 1606  int OrbitalInfoCore(UInt_t run, TFile *f Line 1697  int OrbitalInfoCore(UInt_t run, TFile *f
1697                   bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;                   bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1698               }               }
1699                if(atime>RTpluto2[MU] && atime<=RTstart[MU+1]){                if(atime>RTpluto2[MU] && atime<=RTstart[MU+1]){
1700                   Double_t s_dBdt2=(RTbpluto2[MU] - RTbank2[MU])/(RTpluto2[MU]-RTstart[MU+1]);                   Double_t s_dBdt2=(RTbpluto2[MU] - RTbank2[MU])/(Int_t)(RTpluto2[MU]-RTstart[MU+1]);
1701                   Double_t s_t2=((Double_t)RTpluto2[MU]+(Double_t)RTstart[MU+1])/2. - RTstart[MU];                   Double_t s_t2=((Double_t)RTpluto2[MU]+(Double_t)RTstart[MU+1])/2. - RTstart[MU];
1702                   Double_t s_t1=RTstart[MU+1]-RTstart[MU];                   Double_t s_t1=RTstart[MU+1]-RTstart[MU];
1703                   Double_t s_k=s_dBdt2/(s_t2-s_t1);                   Double_t s_k=s_dBdt2/(s_t2-s_t1);
# Line 1616  int OrbitalInfoCore(UInt_t run, TFile *f Line 1707  int OrbitalInfoCore(UInt_t run, TFile *f
1707                   Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;                   Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1708                 bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;                 bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1709               }               }
1710                orbitalinfo->etha=bank;                if(TMath::Abs(orbitalinfo->etha-bank)>0.1){
1711                Double_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix                  orbitalinfo->etha=bank;
1712                    Double_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1713    
1714                  //Estimations of pitch angle of satellite                  //Estimations of pitch angle of satellite
1715                if(TMath::Abs(bank)>0.7){                  if(TMath::Abs(bank)>0.7){
1716                   Float_t spitch1=TMath::DegToRad()*0.7*diro;//RTdir1[MU];                     Float_t spitch1=TMath::DegToRad()*0.7*diro;//RTdir1[MU];
1717                   Float_t spitch2=TMath::DegToRad()*0.7*diro;//RTdir2[MU];                     Float_t spitch2=TMath::DegToRad()*0.7*diro;//RTdir2[MU];
1718                   Float_t kva=(spitch2-spitch1)/(RTtime2[MU]-RTtime1[MU]);                     Float_t kva=(spitch2-spitch1)/(RTtime2[MU]-RTtime1[MU]);
1719                   Float_t bva=spitch1-kva*RTtime1[MU];                     Float_t bva=spitch1-kva*RTtime1[MU];
1720                   spitch=kva*atime+bva;                     spitch=kva*atime+bva;
1721                }                  }
1722    
1723                //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg                  //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg
1724                Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix                  Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1725                if(TMath::Abs(tlat)<70)                  if(TMath::Abs(tlat)<70)
1726                  yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;                    yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;
1727                yaw = diro*yaw;   //because should be different sign for ascending and descending orbits!                  yaw = diro*yaw; //because should be different sign for ascending and descending orbits!
1728                orbitalinfo->phi=yaw;                  orbitalinfo->phi=yaw;
1729    
1730                if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;                  if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;
1731    
1732  //            Qiji = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch);             // 10RED CHECK  //              Qiji = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch);           // 10RED CHECK
1733                Qij = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch);              // STANDARD                  Qij = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch);            // STANDARD
1734                orbitalinfo->qkind = 1;                  orbitalinfo->qkind = 1;
1735                 }
1736    
1737            //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij            //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij
1738    
1739          } // end of if(atime<RTtime1[0]          } // end of if(atime<RTtime1[0]
1740          } // end of f(((orbitalinfo->TimeGap>60.0 && TMath...          } // end of (((orbitalinfo->TimeGap>60.0 && TMath...
1741         } // end of MU~=0         } // end of MU~=0
1742    
1743          TMatrixD qij = PO->ColPermutation(Qij);          TMatrixD qij = PO->ColPermutation(Qij);
# Line 2121  int OrbitalInfoCore(UInt_t run, TFile *f Line 2214  int OrbitalInfoCore(UInt_t run, TFile *f
2214        //        //
2215        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
2216        //        //
2217          //      tor.Clear("C"); // memory leak?
2218          tor.Delete(); // memory leak?      
2219        delete t_orb;        delete t_orb;
2220        //        //
2221      }; // loop over the events in the run        //      printf(" q0 size %i q0 capacity %i \n",(int)q0.size(),(int)q0.capacity());
2222        } // loop over the events in the run
2223    
2224    
2225      //      //
2226      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
2227      //      //
2228    
2229        //    OrbitalInfotr->FlushBaskets();
2230    
2231      if ( verbose ) printf(" Clear before new run \n");      if ( verbose ) printf(" Clear before new run \n");
2232      delete dbtime;      delete dbtime;
2233    
# Line 2143  int OrbitalInfoCore(UInt_t run, TFile *f Line 2243  int OrbitalInfoCore(UInt_t run, TFile *f
2243      if ( verbose ) printf(" Clear before new run4 \n");      if ( verbose ) printf(" Clear before new run4 \n");
2244      if ( RYPang_lower ) delete RYPang_lower;      if ( RYPang_lower ) delete RYPang_lower;
2245    
2246      if ( l0tr ) l0tr->Delete();  
2247            if ( l0tr ){
2248          if ( verbose ) printf(" delete l0tr\n");
2249          l0tr->Delete();
2250          l0tr = 0;
2251        }
2252        //    if ( l0head ){
2253        //      printf(" delete l0head\n");
2254        //  l0head->Reset();
2255        //  delete l0head;
2256        //  printf(" delete l0head done\n");
2257        //  l0head = 0;
2258        // }
2259        if (eh) {
2260          if ( verbose ) printf(" delete eh\n");
2261          delete eh;
2262          eh = 0;
2263        }
2264    
2265        if ( verbose ) printf(" close file \n");
2266        if ( l0File ) l0File->Close("R");
2267      if ( verbose ) printf(" End run \n");      if ( verbose ) printf(" End run \n");
2268    
2269    }; // process all the runs      q0.clear();
2270          q1.clear();
2271        q2.clear();
2272        q3.clear();
2273        qtime.clear();
2274        qPitch.clear();
2275        qRoll.clear();
2276        qYaw.clear();
2277        qmode.clear();
2278    
2279        if (ch){
2280          if ( verbose ) printf(" delete ch\n");
2281          ch->Delete();
2282          ch = 0;
2283        }
2284      } // process all the runs    <===
2285      if ( debug ){
2286        printf("AFTER LOOP ON RUNs\n");
2287        gObjectTable->Print();
2288      }
2289      //  
2290    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
2291    //    //
2292   closeandexit:   closeandexit:
# Line 2215  int OrbitalInfoCore(UInt_t run, TFile *f Line 2353  int OrbitalInfoCore(UInt_t run, TFile *f
2353    if ( gltle ) delete gltle;    if ( gltle ) delete gltle;
2354    if ( glparam ) delete glparam;    if ( glparam ) delete glparam;
2355    if ( glparam2 ) delete glparam2;    if ( glparam2 ) delete glparam2;
   if ( glparam3 ) delete glparam3;  
2356    if (verbose) printf("\n Exiting3...\n");    if (verbose) printf("\n Exiting3...\n");
2357    if ( glroot ) delete glroot;    if ( glroot ) delete glroot;
2358    if (verbose) printf("\n Exiting4...\n");    if (verbose) printf("\n Exiting4...\n");
# Line 2268  int OrbitalInfoCore(UInt_t run, TFile *f Line 2405  int OrbitalInfoCore(UInt_t run, TFile *f
2405    //    //
2406    if ( debug )  file->ls();    if ( debug )  file->ls();
2407    //    //
2408      if ( debug ){
2409        printf("BEFORE EXITING\n");
2410        gObjectTable->Print();
2411      }
2412    if(code < 0)  throw code;    if(code < 0)  throw code;
2413    return(code);    return(code);
2414  }  }
# Line 2374  void inclresize(vector<Double_t>& t,vect Line 2515  void inclresize(vector<Double_t>& t,vect
2515    
2516  // geomagnetic calculation staff  // geomagnetic calculation staff
2517    
 //void GM_ScanIGRF(TString PATH, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)  
2518  void GM_ScanIGRF(TSQLServer *dbc, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)  void GM_ScanIGRF(TSQLServer *dbc, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
2519  {  {
2520    GL_PARAM *glp = new GL_PARAM();    GL_PARAM *glp = new GL_PARAM();
# Line 2422  void GM_ScanIGRF(TSQLServer *dbc, GMtype Line 2562  void GM_ScanIGRF(TSQLServer *dbc, GMtype
2562          H1->element[23] = temp * 5 + H1->element[22];          H1->element[23] = temp * 5 + H1->element[22];
2563          H1->element[24] = temp * 5 + H1->element[23];          H1->element[24] = temp * 5 + H1->element[23];
2564    if ( glp ) delete glp;    if ( glp ) delete glp;
2565      /*
2566      printf("############################## SCAN IGRF ######################################\n");
2567      printf("       G0      G1     H1\n");
2568      printf(" size  %10i %10i %10i \n",G0->size,G1->size,H1->size);
2569      for ( i = 0; i < 30; i++){
2570        printf("%5i  %10.2f %10.2f %10.2f \n",i,G0->element[i],G1->element[i],H1->element[i]);
2571      }  
2572      printf("###############################################################################\n");
2573      */
2574    } /*GM_ScanIGRF*/
2575    
2576    
2577    
2578    
2579    void GM_SetIGRF(Int_t isSecular, TString ifile1, TString ifile2, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
2580    {
2581      /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/
2582      int i;
2583      double temp,temp2;
2584      int it1,it2;
2585      char buffer[200];
2586      FILE *IGRF;
2587      G0->size = 2;
2588      G1->size = 2;
2589      H1->size = 2;
2590    
2591      for( i = 0; i < 30; i++){
2592        G0->element[i] = 0.;
2593        G1->element[i] = 0.;
2594        H1->element[i] = 0.;
2595      }
2596    
2597      IGRF = fopen(ifile1.Data(), "r");
2598      for( i = 0; i < 2; i++){
2599        fgets(buffer, 200, IGRF);
2600      }
2601      fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G0->element[0],&temp);
2602      fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G1->element[0],&H1->element[0]);
2603      fclose(IGRF);
2604    
2605      IGRF = fopen(ifile2.Data(), "r");
2606      for( i = 0; i < 2; i++){
2607        fgets(buffer, 200, IGRF);
2608      }
2609      if ( isSecular ){
2610        fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2,&temp,&temp2);
2611        G0->element[1] = temp * 5. + G0->element[0];
2612        fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2,&temp,&temp2);
2613        G1->element[1] = temp * 5. + G1->element[0];
2614        H1->element[1] = temp2 * 5. + H1->element[0];
2615      } else {
2616        fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G0->element[1],&temp);
2617        fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G1->element[1],&H1->element[1]);
2618      }
2619      fclose(IGRF);
2620      /*
2621      printf("############################## SCAN IGRF ######################################\n");
2622      printf("       G0      G1     H1\n");
2623      printf(" size  %10i %10i %10i \n",G0->size,G1->size,H1->size);
2624      for ( i = 0; i < 30; i++){
2625        printf("%5i  %10.2f %10.2f %10.2f \n",i,G0->element[i],G1->element[i],H1->element[i]);
2626      }  
2627      printf("###############################################################################\n");
2628      */
2629  } /*GM_ScanIGRF*/  } /*GM_ScanIGRF*/
2630    
2631  void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)  void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)
# Line 2511  void GM_TimeAdjustCoefs(Float_t year, Fl Line 2715  void GM_TimeAdjustCoefs(Float_t year, Fl
2715          int index;          int index;
2716          double x;          double x;
2717          index = (year - GM_STARTYEAR) / 5;          index = (year - GM_STARTYEAR) / 5;
2718          x = (jyear - GM_STARTYEAR) / 5;          x = (jyear - GM_STARTYEAR) / 5.;
2719          Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);          Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);
2720          Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);          Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);
2721          Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);          Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);
# Line 2522  double GM_LinearInterpolation(double x1, Line 2726  double GM_LinearInterpolation(double x1,
2726          /*This function takes a linear interpolation between two given points for x*/          /*This function takes a linear interpolation between two given points for x*/
2727          double weight, y;          double weight, y;
2728          weight  = (x - x1) / (x2 - x1);          weight  = (x - x1) / (x2 - x1);
2729          y = y1 * (1 - weight) + y2 * weight;          y = y1 * (1. - weight) + y2 * weight;
2730            //        printf(" x1 %f x2 %f y1 %f y2 %f x %f ==> y %f \n",x1,x2,y1,y2,x,y);
2731          return y;          return y;
2732  }/*GM_LinearInterpolation*/  }/*GM_LinearInterpolation*/
2733    
Legend:
Removed from v.1.77  
changed lines
  Added in v.1.84
  ViewVC Help
Powered by ViewVC 1.1.23