123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207 |
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- #include "Lemma"
- using namespace Lemma;
-
- #ifdef LEMMAUSEVTK
- #include "matplot.h"
- using namespace matplot;
- #endif
-
- int main() {
-
-
- int nfreq = 6000;
-
-
- int ntimes = 30;
- VectorXr times(ntimes);
- times <<
- 36.0, 45.25, 57.0, 72.25, 92.0,
- 117.0, 148.0, 186.5, 234.0, 290.0,
- 352.5, 427.5, 525.0, 647.5, 802.5,
- 1002.5, 1257.5, 1582.5, 1997.5, 2525.0,
- 3197.5, 4055.0, 5147.5, 6542.5, 8322.5,
- 10592.0, 13490.0, 17187.0, 21902.0, 27915.0;
-
- times = times.array()*1.0e-6;
-
-
-
-
- gaussianquadrature *lgqw = gaussianquadrature::New();
-
- Real lowb;
- Real upb;
- lowb = 1.0;
- upb = 6000000.0;
- lgqw->SetFreqs(nfreq,upb,lowb);
- lgqw->CalcAW();
-
- VectorXr xu(nfreq);
- xu=lgqw->GetAbscissae();
- DipoleSource* Trans=DipoleSource::New();
- Trans->SetType(MAGNETICDIPOLE);
- Trans->SetPolarisation(ZPOLARISATION);
- Trans->SetLocation(0.0,0.0,-1e-4);
-
-
-
-
-
-
-
-
-
-
- Trans->SetNumberOfFrequencies(nfreq);
- VectorXr f(nfreq);
- for (int ifreq=0; ifreq<nfreq; ++ifreq) {
-
- Trans->SetFrequency(ifreq,xu(ifreq));
- f(ifreq) = xu(ifreq);
- }
-
-
-
-
-
- LayeredEarthEM *Earth = LayeredEarthEM::New();
- Earth->SetNumberOfLayers(4);
- Earth->SetLayerConductivity( (VectorXcr(4) << 0.,1.e-6,1.e-2,1.e-6 ).finished() );
- Earth->SetLayerThickness( (VectorXr(2) << 50, 20).finished() );
-
-
- ReceiverPoints *Receivers = ReceiverPoints::New();
- Vector3r loc;
- Real ox = 50.;
- Real oy = 50.;
- Real depth = -1e-2;
- Receivers->SetNumberOfReceivers(1);
- loc << ox, oy, depth;
- Receivers->SetLocation(0, loc);
-
-
- EMEarth1D *EmEarth = EMEarth1D::New();
-
- EmEarth->AttachDipoleSource(Trans);
- EmEarth->AttachLayeredEarthEM(Earth);
- EmEarth->AttachReceiverPoints(Receivers);
- EmEarth->SetFieldsToCalculate(H);
-
- EmEarth->SetHankelTransformMethod(CHAVE);
-
-
-
-
-
- EmEarth->MakeCalc3();
-
-
-
- VectorXcr Hw(nfreq);
- VectorXr e(nfreq);
- Hw(0) = 0.;
- for (int ifreq=0; ifreq<nfreq; ++ifreq) {
- Hw(ifreq) = (Receivers->GetHfield(ifreq, 0))(2);
- e(ifreq) = ifreq;
- }
-
- VectorXr Ht(ntimes);
-
-
- VectorXr func(nfreq);
- for (int ii=0;ii<ntimes;++ii) {
- func = (Hw.imag().array()*((f.array()*times(ii)).sin()));
- lgqw->SetFunc(func);
- lgqw->Integrate();
- Ht(ii)=lgqw->GetResult();
- }
- Ht = Ht.array()*(-2.0/PI);
-
- std::ofstream myfile1;
- myfile1.open("solution.txt");
- for (int ii=0;ii<ntimes;++ii) {
- myfile1 << times(ii) << " " << Ht(ii) << std::endl;
- }
- myfile1.close();
-
- std::ofstream myfile2;
- myfile2.open("freq_domain.txt");
- for (int ii=0;ii<nfreq;++ii) {
- myfile2 << f(ii) << " " << Hw(ii).real() << " " << Hw(ii).imag()
- << std::endl;
- }
- myfile2.close();
-
-
-
-
-
- #ifdef LEMMAUSEVTK
- double colour1[3] = {0.0,0.0,1.0};
- double colour2[3] = {1.0,0.0,0.0};
-
- Plot2D_VTK p1("Hz", "Re(Hw)", 800, 600);
- p1.plot(f, Hw.real().eval(), colour1, ".-");
- p1.show();
-
- Plot2D_VTK p2("Hz", "Im(Hw)", 800, 600);
- p2.plot(f, Hw.imag().eval(), colour2, ".-");
- p2.show();
-
- Plot2D_VTK p3("Time (seconds)", "Ht", 800, 600);
- p3.plot(times, Ht, colour2, ".-");
- p3.show();
-
-
-
-
-
- #endif
-
-
- Trans->Delete();
- Earth->Delete();
- Receivers->Delete();
- EmEarth->Delete();
-
- return EXIT_SUCCESS;
- }
|