/[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.91 by pam-ts, Mon Apr 6 09:22:26 2015 UTC revision 1.92 by pamela, Tue Nov 17 10:10:06 2015 UTC
# Line 1180  int OrbitalInfoCore(UInt_t run, TFile *f Line 1180  int OrbitalInfoCore(UInt_t run, TFile *f
1180        //        //
1181        if ( debug ) printf(" %i start processing \n",procev);              if ( debug ) printf(" %i start processing \n",procev);      
1182        orbitalinfo->Clear();        orbitalinfo->Clear();
1183    
1184        //        //
1185        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
1186        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
# Line 1203  int OrbitalInfoCore(UInt_t run, TFile *f Line 1204  int OrbitalInfoCore(UInt_t run, TFile *f
1204        //        //
1205        if ( debug ) printf(" %i sgp4 \n",procev);              if ( debug ) printf(" %i sgp4 \n",procev);      
1206        cCoordGeo coo;        cCoordGeo coo;
1207        Float_t jyear=0.;            Float_t jyear=0.;
1208        //        //
1209        if(atime >= gltle->GetToTime() || atime < gltle->GetFromTime() ) {  // AGH! bug when reprocessing??        if(atime >= gltle->GetToTime() || atime < gltle->GetFromTime() ) {  // AGH! bug when reprocessing??
1210    
1211          if ( !gltle->Query(atime, dbc) ){          if ( !gltle->Query(atime, dbc) ){
1212            //                  //      
1213            // Compute the magnetic dipole moment.            // Compute the magnetic dipole moment.
# Line 1261  int OrbitalInfoCore(UInt_t run, TFile *f Line 1263  int OrbitalInfoCore(UInt_t run, TFile *f
1263                  for (UInt_t ui = 0; ui < 6; ui++){                  for (UInt_t ui = 0; ui < 6; ui++){
1264                    if (ui>0){                    if (ui>0){
1265                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
                       if ( debug ) printf(" here1 %i \n",ui);  
1266                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
1267                        Int_t recSize = recqtime.size();                        Int_t recSize = recqtime.size();
1268                        if(lowerqtime > recqtime[recSize-1]){                        if(lowerqtime > recqtime[recSize-1]){
# Line 1324  int OrbitalInfoCore(UInt_t run, TFile *f Line 1325  int OrbitalInfoCore(UInt_t run, TFile *f
1325                        }                        }
1326                      }                      }
1327                    }else{                    }else{
1328                      if ( debug ) printf(" here2 %i \n",ui);  //                    if ( debug ) printf(" here2 %i \n",ui);
1329                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
1330                      if(lowerqtime>u_time)nt=0;                      if(lowerqtime>u_time)nt=0;
1331                      Int_t recSize = recqtime.size();                      Int_t recSize = recqtime.size();
# Line 1363  int OrbitalInfoCore(UInt_t run, TFile *f Line 1364  int OrbitalInfoCore(UInt_t run, TFile *f
1364                             qRoll[sizeqmcmd]=RYPang_upper->Kren;                             qRoll[sizeqmcmd]=RYPang_upper->Kren;
1365                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1366                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
 /* CHECK RECOVERED QUATERNIONS PROBLEM */  
 if(recqtime[mu]>=1160987921.75 && recqtime[mu]<=1160987932.00){  
   cout<<"We found it while checking all quaternions"<<"\t"<<recqtime[mu]<<endl;  
 }  
1367                           }                           }
1368                        }                        }
1369                        if(recqtime[mu]>=u_time){                        if(recqtime[mu]>=u_time){
# Line 1394  if(recqtime[mu]>=1160987921.75 && recqti Line 1391  if(recqtime[mu]>=1160987921.75 && recqti
1391                  }                  }
1392                }                }
1393              }              }
             //if ( debug ) cout << "subpacket " << j3 << "\t qtime = " << qtime[qtime.size()-1] << endl;  
1394            }            }
1395          }          }
1396                    
# Line 1420  if(recqtime[mu]>=1160987921.75 && recqti Line 1416  if(recqtime[mu]>=1160987921.75 && recqti
1416          }          }
1417                    
1418    
1419          if ( debug ) printf(" puffi \n");  //      if ( debug ) printf(" puffi \n");
1420          Double_t tmin = 9999999999.;          Double_t tmin = 9999999999.;
1421          Double_t tmax = 0.;          Double_t tmax = 0.;
1422          for(UInt_t tre = 0;tre<qtime.size();tre++){          for(UInt_t tre = 0;tre<qtime.size();tre++){
# Line 1476  if(recqtime[mu]>=1160987921.75 && recqti Line 1472  if(recqtime[mu]>=1160987921.75 && recqti
1472        //Filling Inclination information        //Filling Inclination information
1473        Double_t incli = 0;        Double_t incli = 0;
1474        if ( qtime.size() > 1 ){        if ( qtime.size() > 1 ){
1475          if ( debug ) cout << "ok quaternions is exist and mu = " << must << endl;          if(atime<qtime[0]){
1476          if ( debug ) cout << "qtimes[ " << qtime[0] << " , " << qtime[qtime.size()-1] << " ]\tatime = "<<atime<<endl;            for(UInt_t mu = 1;mu<qtime.size()-1;mu++){
1477                if(qtime[mu]>qtime[0]){
1478                  incli = (qPitch[mu]-qPitch[0])/(qtime[mu]-qtime[0]);
1479                  orbitalinfo->theta =  incli*atime+qPitch[mu]-incli*qtime[mu];
1480                  incli = (qRoll[mu]-qRoll[0])/(qtime[mu]-qtime[0]);
1481                  orbitalinfo->etha =  incli*atime+qRoll[mu]-incli*qtime[mu];
1482                  incli = (qYaw[mu]-qYaw[0])/(qtime[mu]-qtime[0]);
1483                  orbitalinfo->phi =  incli*atime+qYaw[mu]-incli*qtime[mu];
1484                  
1485                  incli = (q0[mu]-q0[0])/(qtime[mu]-qtime[0]);
1486                  orbitalinfo->q0 =  incli*atime+q0[mu]-incli*qtime[mu];
1487                  incli = (q1[mu]-q1[0])/(qtime[mu]-qtime[0]);
1488                  orbitalinfo->q1 =  incli*atime+q1[mu]-incli*qtime[mu];
1489                  incli = (q2[mu]-q2[0])/(qtime[mu]-qtime[0]);
1490                  orbitalinfo->q2 =  incli*atime+q2[mu]-incli*qtime[mu];
1491                  incli = (q3[mu]-q3[0])/(qtime[mu]-qtime[0]);
1492                  orbitalinfo->q3 =  incli*atime+q3[mu]-incli*qtime[mu];
1493                  orbitalinfo->TimeGap=qtime[0]-atime;
1494                  break;
1495                }
1496              }
1497            }
1498            Float_t eend=qtime.size()-1;
1499            if(atime>qtime[eend]){
1500              for(UInt_t mu=eend-1;mu>=0;mu--){
1501                if(qtime[mu]<qtime[eend]){
1502                  incli = (qPitch[eend]-qPitch[mu])/(qtime[eend]-qtime[mu]);
1503                  orbitalinfo->theta =  incli*atime+qPitch[eend]-incli*qtime[eend];
1504                  incli = (qRoll[eend]-qRoll[mu])/(qtime[eend]-qtime[mu]);
1505                  orbitalinfo->etha =  incli*atime+qRoll[eend]-incli*qtime[eend];
1506                  incli = (qYaw[eend]-qYaw[mu])/(qtime[eend]-qtime[mu]);
1507                  orbitalinfo->phi =  incli*atime+qYaw[eend]-incli*qtime[eend];
1508              
1509                  incli = (q0[eend]-q0[mu])/(qtime[eend]-qtime[mu]);
1510                  orbitalinfo->q0 =  incli*atime+q0[eend]-incli*qtime[eend];
1511                  incli = (q1[eend]-q1[mu])/(qtime[eend]-qtime[mu]);
1512                  orbitalinfo->q1 =  incli*atime+q1[eend]-incli*qtime[eend];
1513                  incli = (q2[eend]-q2[mu])/(qtime[eend]-qtime[mu]);
1514                  orbitalinfo->q2 =  incli*atime+q2[eend]-incli*qtime[eend];
1515                  incli = (q3[eend]-q3[mu])/(qtime[eend]-qtime[mu]);
1516                  orbitalinfo->q3 =  incli*atime+q3[eend]-incli*qtime[eend];
1517                  break;
1518                }
1519              }
1520            }
1521          for(UInt_t mu = must;mu<qtime.size()-1;mu++){          for(UInt_t mu = must;mu<qtime.size()-1;mu++){
1522            if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);            if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);
1523            if(qtime[mu+1]>qtime[mu]){            if(qtime[mu+1]>qtime[mu]){
# Line 1502  if(recqtime[mu]>=1160987921.75 && recqti Line 1542  if(recqtime[mu]>=1160987921.75 && recqti
1542                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];
1543                Float_t tg = (qtime[mu+1]-qtime[mu])/1000.0;                Float_t tg = (qtime[mu+1]-qtime[mu])/1000.0;
1544                if(tg>=1) tg=0.00;                if(tg>=1) tg=0.00;
1545                orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1])-atime,TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];                orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1]-atime),TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];
1546                orbitalinfo->mode = qmode[mu+1];                orbitalinfo->mode = qmode[mu+1];
1547                //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;                //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;
1548                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;
# Line 1518  if(recqtime[mu]>=1160987921.75 && recqti Line 1558  if(recqtime[mu]>=1160987921.75 && recqti
1558        //        //
1559                
1560        if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){        if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){
1561          if ( debug ) cout << "ops no iclination information" << endl;          if (debug) cout << "Warning: no iclination information "<< endl;
1562          orbitalinfo->mode = 10;          orbitalinfo->mode = 10;
1563          orbitalinfo->q0 = -1000.;          orbitalinfo->q0 = -1000.;
1564          orbitalinfo->q1 = -1000.;          orbitalinfo->q1 = -1000.;
# Line 1528  if(recqtime[mu]>=1160987921.75 && recqti Line 1568  if(recqtime[mu]>=1160987921.75 && recqti
1568          orbitalinfo->phi = -1000.;          orbitalinfo->phi = -1000.;
1569          orbitalinfo->theta = -1000.;          orbitalinfo->theta = -1000.;
1570          orbitalinfo->TimeGap = -1000.;          orbitalinfo->TimeGap = -1000.;
1571            TMatrixD Iij(3,3);
1572            Iij(0,0)=0; Iij(0,1)=0; Iij(0,2)=0;
1573            Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
1574            Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
1575            orbitalinfo->Iij.ResizeTo(Iij);
1576            orbitalinfo->Iij = Iij;
1577            
1578          //orbitalinfo->qkind = -1000;          //orbitalinfo->qkind = -1000;
1579                    
1580          //      if ( debug ){          //      if ( debug ){
# Line 1586  if(recqtime[mu]>=1160987921.75 && recqti Line 1633  if(recqtime[mu]>=1160987921.75 && recqti
1633                
1634        if ( debug ) printf(" coord done \n");        if ( debug ) printf(" coord done \n");
1635        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
         //        
1636          orbitalinfo->lon = lon;          orbitalinfo->lon = lon;
1637          orbitalinfo->lat = lat;          orbitalinfo->lat = lat;
1638          orbitalinfo->alt = alt;          orbitalinfo->alt = alt;
# Line 1646  if(recqtime[mu]>=1160987921.75 && recqti Line 1692  if(recqtime[mu]>=1160987921.75 && recqti
1692          //14.9/(xl*xl);          //14.9/(xl*xl);
1693          orbitalinfo->igrf_icode = (Float_t)icode;          orbitalinfo->igrf_icode = (Float_t)icode;
1694          //          //
1695        }              }  //check lon lat alt      
1696        //        //
1697        if ( debug ) printf(" pitch angle \n");        if ( debug ) printf(" pitch angle \n");
1698        //        //
1699        // pitch angles        // pitch angles
1700        //        //
1701        if( orbitalinfo->TimeGap>0){        if( orbitalinfo->TimeGap>=0){
1702          //          //
1703          if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);          if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);
1704          Float_t Bx = -orbitalinfo->Bdown;          Float_t Bx = -orbitalinfo->Bdown;
# Line 1673  if(recqtime[mu]>=1160987921.75 && recqti Line 1719  if(recqtime[mu]>=1160987921.75 && recqti
1719          if (orbitalinfo->TimeGap>180.0) tgpar0=true;          if (orbitalinfo->TimeGap>180.0) tgpar0=true;
1720         if(MU!=0){         if(MU!=0){
1721  //      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)
1722    
1723         if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && (TMath::Abs(orbitalinfo->etha)>0.1 || tgpar0)) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || tgpar))){         if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && (TMath::Abs(orbitalinfo->etha)>0.1 || tgpar0)) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || tgpar))){
1724    
1725          //found in Rotation Table this data for this time interval          //found in Rotation Table this data for this time interval
1726          if(atime<RTtime1[0])          if(atime<RTtime1[0])
1727            orbitalinfo->azim = 5;        //means that RotationTable no started yet            orbitalinfo->azim = 5;        //means that RotationTable no started yet
1728         else{         else{
1729                  // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position                  // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position
1730    
1731                Double_t bank=RTstart[MU];                Double_t bank=RTstart[MU];
1732                Double_t tlat=orbitalinfo->lat;                Double_t tlat=orbitalinfo->lat;
1733    
# Line 1749  if(recqtime[mu]>=1160987921.75 && recqti Line 1798  if(recqtime[mu]>=1160987921.75 && recqti
1798          TMatrixD Gij = PO->ColPermutation(Fij);          TMatrixD Gij = PO->ColPermutation(Fij);
1799          Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);          Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1800          TMatrixD Iij = PO->ColPermutation(Dij);          TMatrixD Iij = PO->ColPermutation(Dij);
1801    
1802          TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);          TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);
1803          // go to Pamela reference frame from Resurs reference frame          // go to Pamela reference frame from Resurs reference frame
1804          Float_t tmpy = SP.Y();          Float_t tmpy = SP.Y();
Legend:
Removed from v.1.91  
changed lines
  Added in v.1.92
  ViewVC Help
Powered by ViewVC 1.1.23