diffangle.cpp

Go to the documentation of this file.

#include "MTStationList.h"
#include "CFatalException.h"
#include <iostream>
#include <string>
#include <fstream>
#include <iomanip>
#include "Jacknife.h"

using namespace std;
using namespace gplib;

string version = "$Id: diffangle.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
void PromptAndReadStationFile(MTStationList &List, std::string prompt)
  {
    string listname;
    cout << prompt;
    cin >> listname;
    try
      {
        List.GetData(listname); //try to read in all station data
      } catch (FatalException &e)
      {
        cerr << e.what() << endl; // if something fails print error
        exit(100); // and stop execution
      }
  }

int PromptAndReadFreqIndex(MTStationList &List, std::string prompt)
  {
    cout << prompt << endl << endl;
    cout << "Available Periods: " << List.GetCommonFrequencies().size() << endl;

    for (trealdata::const_iterator frequencies =
        List.GetCommonFrequencies().begin(); frequencies
        != List.GetCommonFrequencies().end(); ++frequencies)
      {
        cout << setw(9) << setfill(' ') << setprecision(4) << distance(
            List.GetCommonFrequencies().begin(), frequencies);
        cout << setw(15) << setfill(' ') << setprecision(8) << 1.
            / (*frequencies) << endl;
      }
    int freqindex;
    cout << "Select frequency Index: ";
    cin >> freqindex;
    return freqindex;
  }

int main()
  {
    string list1, list2; //storing names for input and root of output
    cout
        << " This is diffangle, calculate the difference in strike angle and phase split"
        << endl; // write some info
    cout << " for corresponding sites in two different lists. " << endl << endl;

    MTStationList MTSites1; // object for the first station list
    PromptAndReadStationFile(MTSites1, "Station List 1:");
    MTStationList MTSites2; // object for the second station list
    PromptAndReadStationFile(MTSites2, "Station List 2:");
    tStatSyncPair StatSync = FindCorrespondingSites(MTSites1.GetList(),
        MTSites2.GetList());
    const std::vector<tindexvector> freqindices1(MTSites1.GetComFreqIndices());
    const std::vector<tindexvector> freqindices2(MTSites2.GetComFreqIndices());
    int SelectedFreq1 = PromptAndReadFreqIndex(MTSites1, "List 1");
    int SelectedFreq2 = PromptAndReadFreqIndex(MTSites2, "List 2");
    ofstream outfile("out");
    const int nbootsamples = 10000;
    CJacknife ErrEst(nbootsamples); //init error estimation object
    for (size_t i = 0; i < StatSync.size(); ++i)
      {

        double strikediff = (StatSync.at(i).first->GetMTData().at(
            freqindices1.at(i).at(SelectedFreq1)).GetPhiEllip()
            - StatSync.at(i).second->GetMTData().at(freqindices2.at(i).at(
                SelectedFreq2)).GetPhiEllip());
        double strikeerr, dummy;
        ErrEst.CalcErrors(&MTTensor::GetPhiEllip,
            StatSync.at(i).first->GetMTData().at(freqindices1.at(i).at(
                SelectedFreq1)), dummy, strikeerr);
        double phidiffdiff = (StatSync.at(i).first->GetMTData().at(
            freqindices1.at(i).at(SelectedFreq1)).GetBeta_phi()
            - StatSync.at(i).second->GetMTData().at(freqindices2.at(i).at(
                SelectedFreq2)).GetBeta_phi());
        double phidifferr;
        ErrEst.CalcErrors(&MTTensor::GetBeta_phi,
            StatSync.at(i).first->GetMTData().at(freqindices1.at(i).at(
                SelectedFreq1)), dummy, phidifferr);
        outfile << phidiffdiff * 180. / PI << " " << strikediff << " "
            << phidifferr * 180. / PI << " " << strikeerr << endl;
      }
  }