00001 #include "AbsVelRecObjective.h"
00002 #include <gsl/gsl_math.h>
00003
00004 AbsVelRecObjective::~AbsVelRecObjective()
00005 {
00006 }
00007
00008
00009 AbsVelRecObjective::AbsVelRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData,
00010 SurfaceWaveData &AbsVel, const int myshift, const double myomega,
00011 const double mysigma, const double myc, const double myslowness,
00012 const RecCalc::trfmethod method) :
00013 C1DRecObjective(TheRecData, myshift, myomega, mysigma, myc, myslowness,
00014 method), CalcAbsVel(myomega,mysigma,myc, method), MeasuredAbsVel(AbsVel)
00015 {
00016 recweight = 1;
00017 absvelweight = 1;
00018 }
00019
00020 AbsVelRecObjective::AbsVelRecObjective(const AbsVelRecObjective &Old) :
00021 C1DRecObjective(Old), CalcAbsVel(Old.CalcAbsVel),
00022 MeasuredAbsVel(Old.MeasuredAbsVel), SynthAbsVel(Old.SynthAbsVel)
00023 {
00024 recweight = Old.recweight;
00025 absvelweight = Old.absvelweight;
00026 }
00027
00028 AbsVelRecObjective& AbsVelRecObjective::operator=(const AbsVelRecObjective& source)
00029 {
00030 if (this == &source) return *this;
00031 C1DRecObjective::operator=(source);
00032 CalcAbsVel = source.CalcAbsVel;
00033 MeasuredAbsVel = source.MeasuredAbsVel;
00034 SynthAbsVel = source.SynthAbsVel;
00035 recweight = source.recweight;
00036 absvelweight = source.absvelweight;
00037 return *this;
00038 }
00039
00040 void AbsVelRecObjective::SafeParallel(const ttranscribed &member)
00041 {
00042 C1DRecObjective::SafeParallel(member);
00043 double returnvalue = pow(GetRMS(), GetFitExponent()) * GetMisfit().size();
00044 returnvalue *= recweight;
00045 if (absvelweight > 0.0 && !MeasuredAbsVel.GetPhaseVelocities().empty())
00046 {
00047 SynthAbsVel.SetPeriods() = MeasuredAbsVel.GetPeriods();
00048 SeismicDataComp Rad(RecCalculator.GetRadComp()),
00049 Ver(RecCalculator.GetVerComp());
00050 CalcAbsVel.CalcRFVel(slowness, Rad, Ver,
00051 SynthAbsVel.SetPhaseVelocities());
00052
00053 assert(SynthAbsVel.GetPhaseVelocities().size() == MeasuredAbsVel.GetPhaseVelocities().size());
00054 const unsigned int nvel = SynthAbsVel.GetPhaseVelocities().size();
00055 const unsigned int nrec =GetMisfit().size();
00056
00057
00058 tmisfit LocalMisfit(nrec+nvel);
00059 tdata LocalData(nrec+nvel);
00060 copy(GetMisfit().begin(), GetMisfit().end(), LocalMisfit.begin());
00061 copy(GetSynthData().begin(), GetSynthData().end(), LocalData.begin());
00062
00063 double currvalue = 0.0;
00064 for (unsigned int i = 0; i < nvel; ++i)
00065 {
00066 currvalue = SynthAbsVel.GetPhaseVelocities().at(i) - MeasuredAbsVel.GetPhaseVelocities().at(i);
00067 currvalue /= MeasuredAbsVel.GetPhaseVelocities().at(i) * errorlevel;
00068 currvalue = gsl_pow_int(currvalue, GetFitExponent());
00069 LocalMisfit(i+nrec)= fabs(currvalue);
00070 LocalData(i+nrec)= SynthAbsVel.GetPhaseVelocities().at(i);
00071 assert(gsl_finite(currvalue));
00072 returnvalue += absvelweight * fabs(currvalue);
00073 }
00074 SetMisfit(LocalMisfit);
00075 SetSynthData(LocalData);
00076 unsigned int ndata;
00077 double divisor = 1.0;
00078 if (recweight > 0)
00079 {
00080 ndata = nrec+nvel;
00081 divisor = nrec * recweight + nvel*absvelweight;
00082 }
00083 else
00084 {
00085 ndata = nvel;
00086 divisor = ndata * absvelweight;
00087 }
00088 SetRMS(pow(returnvalue/divisor, 1.0/GetFitExponent()));
00089 }
00090
00091 }