VISIBILITY CALCUTALTION

The following program is a simple simulator that calculates orbit, attitude and geolocation geometric properties based on a simple instrument (push broom sensor at a fixed azimuth and a certain numbers of pixels). Observation are above a certain area:

// ------------------------------------------------------------------------
//
//                           C++ Example Program.
//                           FULLY-FEATURED PROGRAM
//
// ------------------------------------------------------------------------

#include <string.h>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>

#include "DataHandlingData.h"
#include "AttitudeDefinitionData.h"
#include "TimeCorrelation.h"
#include "ModelId.h"
#include "OrbitId.h"
#include "OrbitExtra.h"
#include "PointingFunc.h"
#include "SatNomTransId.h"
#include "SatTransId.h"
#include "InstrTransId.h"
#include "Attitude.h"
#include "DemId.h"
#include "AtmosId.h"
#include "Target.h"
#include "Swath.h"
#include "SwathStar.h"
#include "VisibilityData.h"
#include "VisibilityList.h"
#include "VisibilitySegment.h"
#include "SatId.h"
#include "StfFile.h"
#include "SdfFile.h"
#include "ZoneFile.h"
#include "StationFile.h"
#include "ZoneCoverageIn.h"
#include "ZoneCoverageOut.h"

using namespace std;
using namespace EECFI;

// ------------ 
// Main program 
// ------------ 
int main (int argc, char *argv[])
{
 
  vector<string> errorMessages;

  // Variables to initialize dem_id
  long useDEMFlag = 0; // If 0, DEM will NOT be used. If 1, DEM will be used

  try
  {
    // Initialize default model id
    long mode = XLCFI_MODEL_DEFAULT;
    vector<long>models; //not used when mode is XLCFI_MODEL_DEFAULT
    ModelId *modelId = new ModelId(mode, models);
  
    // Initialize time correlations with orbit file 
    // (using all the file validity interval)
    string orbitFile = "orbit_file.xml";
    long timeModel = XLCFI_TIMEMOD_AUTO;
    long timeRef   = XLCFI_TIME_UTC;
    vector <string> inputFiles;
    inputFiles.push_back(orbitFile);
    TimeInterval validityTime;
    
    TimeCorrelation timeId(timeModel, inputFiles, timeRef, validityTime);
   
    // Initialize orbit 
    // The same orbit files as for TimeCorrelation will be used
    SatId satId(XOCFI_SAT_METOP_SG_A1); // METOP satellite
    long   initModel = XOCFI_ORBIT_INIT_AUTO; // Detect automatically the type of orbit file
    long   initMode  = XOCFI_SEL_FILE; // Use all file, time and orbit intervals are dummy
    // time0, time1, orbi0, orbit1 are not used for XLCFI_SEL_FILE option
    double time0     = 0.;
    double time1     = 0.;
    long   orbit0    = 0;
    long   orbit1    = 0;
    
    OrbitId orbitId(satId, *modelId, timeId, 
                    timeRef, initModel, inputFiles, 
                    initMode, time0, time1, orbit0, orbit1);

    // Initialize dem id if required
    DemId demId;
    if (useDEMFlag)
    {
      string demConfigFile = "dem_config_file.xml";
      initMode  = XDCFI_DEM_GETASSE30_V1; // it is not needed, but it is better to have it
                                          // in line with the content of DEM configuration file
      initModel = 0; // Dummy parameter
      
      demId.init(initMode, initModel, demConfigFile);
      
    }
    
    // Initialize atmos id with no atmosphere model
    initMode  = XPCFI_NO_REF_INIT; // No atmosphere model
    initModel = 0;              // dummy parameter
    string atmosDummyFile = ""; // dummy parameter
    AtmosId atmosId;
    atmosId.init(initMode, initModel, atmosDummyFile);

    // Initialize satellite attitudes to perform.
    // We use an attitude definition file, this way we can initialize all the
    // atttitudes at the same time.
    AttitudeDefinitionData attDefData("attitude_definition_file.xml");
    attDefData.read();
    
    AttitudeDef attDef;
    PointingFunc::attitudeDefine(attDefData,
                                 attDef.satNomTransId,
                                 attDef.satTransId,
                                 attDef.instrTransId); 
    attDef.type = XPCFI_SAT_ATT;
    
    Attitude attId(timeId, 
                   attDef.satNomTransId,
                   attDef.satTransId,
                   attDef.instrTransId);


    // Define zone (Spain polygon)
    vector<ZonePoint> zonePoints;
    zonePoints.push_back(ZonePoint( 3.813258057950049,  42.5091355830345));
    zonePoints.push_back(ZonePoint(-1.706802028054565,  35.82814264689165));
    zonePoints.push_back(ZonePoint(-9.155852980327245,  36.05894971586847));
    zonePoints.push_back(ZonePoint(-9.876611106388474,  44.21197011498435));
    zonePoints.push_back(ZonePoint( 3.813258057950049,  42.5091355830345));
    ZoneRec zoneData("TEST_ZONE",
                     "Spain zone", 
                     "", //surface type
                     "User",
                     XDCFI_POLYGON, 
                     XDCFI_RECTANGULAR,
                     0, //diameter
                     zonePoints);
    
    vector<ZoneInfo> zoneVec;
    zoneVec.push_back(ZoneInfo(zoneData, 
                               XDCFI_RECTANGULAR, 
                               0.)); // min duration
    ZoneInfoList zoneInfoList(XVCFI_COMPUTE, zoneVec);

    // Define visibility interval from start of orbit validity and 1 day more.
    ValidityTime orbitVal =     orbitId.getValTime();
    Time valTime0 (timeId, orbitVal.start, orbitVal.timeRef);

    Time valUTC0 = valTime0.getProcessing(XLCFI_TIME_UTC);
    
    VisTimeInterval interval;
    interval.tstart.type    = XVCFI_UTC_TYPE;
    interval.tstart.utcTime = valUTC0.time;
    interval.tstop.type     = XVCFI_UTC_TYPE;
    interval.tstop.utcTime  = valUTC0.time + 1.;

    // Initialize swath id with Swath definition file
    SwathInfo swathInfo;
    swathInfo.type           = XVCFI_FILE_SDF; // Initialize with Swath definition file
    swathInfo.filename       = "swath_definition_file.xml";
    swathInfo.nofRegenOrbits = 10; // Regenerate internal STF every 10 orbits
    //swathInfo.sdfFile = NULL;
    //swathInfo.stfFile = NULL;

    SwathId swathId(atmosId, swathInfo);
    
    Swath swath(orbitId, 1, "");
    VisibilityList visList = swath.zoneVisTimeCompute(attDef,
                                                      swathId,
                                                      zoneInfoList,
                                                      interval);

    cout << "Number of visibility segments found: "<< visList.size() << endl;

    // Prepare inputs for xp_target_list_inter
    long nofPixels = 3;
    long az_i;
    
    InstrumentData instData;
    instData.type            = XPCFI_AZ_EL_LIST;
    instData.signalFrequency = 1e10;

    double azimuthList[3]   = {90.,90.,270.};
    double elevationList[3] = {75.,90.,60.};
    
    for (az_i = 0 ; az_i < nofPixels ; az_i ++)
    {
      AzimuthElevation azEl;
      azEl.azimuth        = azimuthList[az_i];
      azEl.elevation      = elevationList[az_i];
      azEl.azimuthRate   = 0.;
      azEl.elevationRate = 0.;
      instData.azimuthElevationList.azElList.push_back(azEl);
    }
    
  
    // Perform some computations in the visibility intervals
    // We will use a time step of 5 seconds:
    double timeStep = 5. / 86400.;
    
    // Select the type of target to be computed
    long targetType;
    if (useDEMFlag) targetType = XPCFI_DEM_TARGET_TYPE;
    else            targetType = XPCFI_USER_TARGET_TYPE;
    
    // Open csv info file to store information and write header line
    ofstream csvFile;
    csvFile.open ("out.csv");
    csvFile << "UTC_TIME,SAT_POSITION_X,SAT_POSITION_Y,SAT_POSITION_Z," 
            << "SAT_VELOCITY_X,SAT_VELOCITY_Y,SAT_VELOCITY_Z,"
            << "PIXEL_0_TARGET_LON,PIXEL_0_TARGET_LAT,"
            << "PIXEL_1_TARGET_LON,PIXEL_1_TARGET_LAT,"
            << "PIXEL_2_TARGET_LON,PIXEL_2_TARGET_LAT"
            << "\n";

    long seg_i;
    long propagModel = 0; // Dummy

    for (seg_i = 0 ; seg_i < visList.size() ; seg_i ++)
    {
      // Loop to visibility interval
      double timePointsInInterval = (visList[seg_i].utcTimeStop - visList[seg_i].utcTimeStart) / timeStep;
      
      long time_i;
      for (time_i = 0 ; time_i < timePointsInInterval; time_i ++)
      {
        double loopTime = visList[seg_i].utcTimeStart + time_i * timeStep;
                    
        // Compute Orbit State Vector at current time
        StateVector osv =       orbitId.osvCompute(propagModel, timeRef, loopTime);
            
        // Compute attitude at current time (satellite attitude)
        attId.compute (*modelId, osv, XPCFI_SAT_ATT);
        
        // Compute target METHOD 1: Compute pixels one by one
        long pixel_i;
        for (pixel_i = 0 ; pixel_i < nofPixels ; pixel_i++)
        {
          Target targetId(attId, demId, atmosId);
          
          targetId.targetInter(satId, 
                               XPCFI_NO_DER, 
                               XPCFI_INTER_1ST,
                               azimuthList[pixel_i], 
                               elevationList[pixel_i],
                               0.,     //geodetic altitude 
                               0., 0,  //azimuth and elevation rates
                               0,      //iray (dummy)
                               1.e10); //signal frequency

          vector<double> mainResult;
          vector<double> mainResultRate;
          vector<double> mainResultRateRate;

          targetId.extraMain(XPCFI_TARG_EXTRA_MAIN_GEO, 
                             XPCFI_USER_TARGET_TYPE, 
                             0, // targetNumber, 
                             mainResult,
                             mainResultRate,
                             mainResultRateRate);
                             
        } // end loop to pixels
        
        
        // Compute target METHOD 2: Compute all pixels at the same time
        Target targetListId(attId, demId, atmosId);
        TargetOutput targOut;
        targetListId.targetListInter(satId,
                                     XPCFI_NO_DER, 
                                     XPCFI_INTER_1ST, 
                                     instData, 
                                     0.,    //geodetic altitude
                                     targOut);

        TargetExtraResultsList resultslist;                                     
        targetListId.extraListMain(XPCFI_TARG_EXTRA_MAIN_GEO, 
                                   XPCFI_USER_TARGET_TYPE,
                                   resultslist);
        
        csvFile << loopTime << ","
                << osv.coord.pos[0] << "," << osv.coord.pos[1] << "," << osv.coord.pos[2] << ","
                << osv.coord.vel[0] << "," << osv.coord.vel[1] << "," << osv.coord.vel[2] << ","
                << resultslist.extraResults[0].results[0] << "," 
                << resultslist.extraResults[0].results[2] << ","
                << resultslist.extraResults[1].results[0] << "," 
                << resultslist.extraResults[1].results[2] << ","
                << resultslist.extraResults[2].results[0] << "," 
                << resultslist.extraResults[2].results[2] << ","
                << endl;
                
      } // End loop to interval
      
    } // end loop to visibility intervals

    // Close output file
    csvFile.close();
    
  }
  catch( CfiError cfiError )
  {
    cout << "\tERROR IN EXAMPLE!!!" << endl;
    errorMessages = cfiError.getMsg( errorMessages );
  }

  //Print result
  if ( errorMessages.size() > 0)
  {
    cout << endl << endl << "*** ERRORS DURING EXECUTION ***" << endl;
    for (int cont = 0 ; cont < errorMessages.size() ; cont++ )
    {
      cout << endl << "Message error " << cont + 1 << ":"
      << errorMessages[cont] << endl;
    }
    cout << "*** EXAMPLE: FAILED ***" << endl << endl;
    
  }//end if errVec.size
  else
  {
    cout << endl << "*** USAGE EXAMPLE COMPLETED SUCCESSFULLY ***" << endl << endl;
  }
  
  return(0);
}