Browse Source

Key201 working, but needs cleanup

lagkey
Trevor Irons 6 years ago
parent
commit
1eebe3dfc3

+ 1
- 0
Modules/FDEM1D/include/HankelTransform.h View File

42
 
42
 
43
                 /// @todo this is actually a deprecated function. Only Chave
43
                 /// @todo this is actually a deprecated function. Only Chave
44
                 ///   uses this.
44
                 ///   uses this.
45
+                /// @deprecated
45
                 /// Performs numerical integration using Gaussian quadrature
46
                 /// Performs numerical integration using Gaussian quadrature
46
                 /// ikk:   type of kernel depending on source and receiver couple
47
                 /// ikk:   type of kernel depending on source and receiver couple
47
                 /// imode: a switch for TE(0) and TM(1) mode
48
                 /// imode: a switch for TE(0) and TM(1) mode

+ 9
- 4
Modules/FDEM1D/src/EMEarth1D.cpp View File

778
         }
778
         }
779
 
779
 
780
         // Determine number of lagged convolutions to do
780
         // Determine number of lagged convolutions to do
781
-        // TODO, can Hankel2 adjust the lagg spacing safely?
782
-        int nlag = 1; // We need an extra for some reason for stability
781
+        int nlag = 1; // (Key==0)  We need an extra for some reason for stability? Maybe in Spline?
783
         Real lrho ( 1.01* rhomax );
782
         Real lrho ( 1.01* rhomax );
784
         while ( lrho > rhomin ) {
783
         while ( lrho > rhomin ) {
785
             nlag += 1;
784
             nlag += 1;
786
             lrho *= Hankel->GetABSER();
785
             lrho *= Hankel->GetABSER();
786
+            //std::cout << "lrho\t" << lrho << std::endl;
787
         }
787
         }
788
+        // Key variant
789
+//        Real lamMin = Filter.base(1)/rMax;
790
+//        Real lamMax = Filter.base(end)/rMin;
791
+//        auto nLambda = ceil(log(lamMax/lamMin)/log(filterSpacing))+1;
792
+        //nlag = 3;
793
+        //std::cout << "nlag\t" << nlag << std::endl;
788
 
794
 
789
         //int nlag = rhomin
795
         //int nlag = rhomin
790
         auto tDipole = antenna->GetDipoleSource(0);
796
         auto tDipole = antenna->GetDipoleSource(0);
793
         // Instead we should pass the antenna into this so that Hankel hass all the rho arguments...
799
         // Instead we should pass the antenna into this so that Hankel hass all the rho arguments...
794
         Hankel->ComputeLaggedRelated( 1.01* rhomax, nlag, tDipole->GetKernelManager() );
800
         Hankel->ComputeLaggedRelated( 1.01* rhomax, nlag, tDipole->GetKernelManager() );
795
 
801
 
796
-        //std::cout << Hankel->GetAnswer() << std::endl;
797
-        //std::cout << Hankel->GetArg() << std::endl;
802
+        //std::cout << "After! " << Hankel->GetAnswer() << std::endl;
798
 
803
 
799
         // Sort the dipoles by rho
804
         // Sort the dipoles by rho
800
         for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
805
         for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {

+ 2
- 2
Modules/FDEM1D/src/FHTAnderson801.cpp View File

950
 
950
 
951
         // Now do cubic spline
951
         // Now do cubic spline
952
         // TODO Check that knots are set in right order, Eigen has reverse()
952
         // TODO Check that knots are set in right order, Eigen has reverse()
953
-        //std::cout << "Arg\n" << Arg << std::endl;
954
-        //std::cout << "Zans\n" << Zans.col(0) << std::endl;
953
+        std::cout << "Arg\n" << Arg << std::endl;
954
+        std::cout << "Zans\n" << Zans << std::endl;
955
         for (int ii=0; ii<Zans.cols(); ++ii) {
955
         for (int ii=0; ii<Zans.cols(); ++ii) {
956
             auto Spline = CubicSplineInterpolator::NewSP();
956
             auto Spline = CubicSplineInterpolator::NewSP();
957
             Spline->SetKnots( Arg, Zans.col(ii).real() );
957
             Spline->SetKnots( Arg, Zans.col(ii).real() );

+ 38
- 12
Modules/FDEM1D/src/FHTKey201.cpp View File

367
         for (int ir2=0; ir2<nrel; ++ir2) {
367
         for (int ir2=0; ir2<nrel; ++ir2) {
368
             Zans(0, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
368
             Zans(0, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
369
         }
369
         }
370
-
370
+        std::cout << "rho\n" << rho << std::endl;
371
+        std::cout << "Zans\n" << Zans << std::endl;
372
+        exit(EXIT_SUCCESS);
371
         return ;
373
         return ;
372
     }		// -----  end of method FHTKey201::ComputeRelated  -----
374
     }		// -----  end of method FHTKey201::ComputeRelated  -----
373
 
375
 
376
     //      Method:  ComputeLaggedRelated
378
     //      Method:  ComputeLaggedRelated
377
     //--------------------------------------------------------------------------------------
379
     //--------------------------------------------------------------------------------------
378
     void FHTKey201::ComputeLaggedRelated ( const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager ) {
380
     void FHTKey201::ComputeLaggedRelated ( const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager ) {
381
+
382
+        //std::cout << "rho max\t " << rho << std::endl;
383
+        //Real rho = 214.963;
384
+
379
         int nrel = (int)(KernelManager->GetSTLVector().size());
385
         int nrel = (int)(KernelManager->GetSTLVector().size());
380
-        Eigen::Matrix<Complex, 201, Eigen::Dynamic > Zwork;
386
+
387
+        Eigen::Matrix< Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
388
+        //Eigen::Matrix<Complex, 201+nrel, Eigen::Dynamic > Zwork;
381
         Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nlag, nrel);
389
         Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nlag, nrel);
382
-        Zwork.resize(201, nrel);
383
-        VectorXr lambda = WT201.col(0).array()/rho;
390
+        Zwork.resize(201+nlag, nrel);  // Zwork needs to be expanded to filter length + nlag
391
+
392
+        // lambda needs to be expanded to include lagged results
393
+        VectorXr lambda = (VectorXr(201+nlag) << WT201.col(0).array()/rho, VectorXr::Zero(nlag)).finished();
394
+        for (int ilam =201; ilam< nlag+201; ++ilam) {
395
+            lambda(ilam) = lambda(ilam-1)/GetABSER();
396
+        }
397
+
384
         int NumFun = 0;
398
         int NumFun = 0;
385
         int idx = 0;
399
         int idx = 0;
386
 
400
 
401
+        //std::cout << lambda.transpose() << std::endl;
402
+
387
         VectorXr Arg(nlag);
403
         VectorXr Arg(nlag);
388
         Arg(nlag-1) = rho;
404
         Arg(nlag-1) = rho;
389
         for (int ilag=nlag-2; ilag>=0; --ilag) {
405
         for (int ilag=nlag-2; ilag>=0; --ilag) {
390
             Arg(ilag) = Arg(ilag+1) * GetABSER();
406
             Arg(ilag) = Arg(ilag+1) * GetABSER();
391
         }
407
         }
408
+        //std::cout << "Arg\t" << Arg << std::endl;
392
 
409
 
393
         // Get Kernel values
410
         // Get Kernel values
394
         for (int ir=0; ir<lambda.size(); ++ir) {
411
         for (int ir=0; ir<lambda.size(); ++ir) {
397
             KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
414
             KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
398
             for (int ir2=0; ir2<nrel; ++ir2) {
415
             for (int ir2=0; ir2<nrel; ++ir2) {
399
                 // Zwork* needed due to sign convention of filter weights
416
                 // Zwork* needed due to sign convention of filter weights
417
+ 			    //Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
400
  			    Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
418
  			    Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
401
             }
419
             }
402
-
403
         }
420
         }
404
 
421
 
405
         // We diverge slightly from Key here, each kernel is evaluated seperately, whereby instead
422
         // We diverge slightly from Key here, each kernel is evaluated seperately, whereby instead
407
         // in the interests of making them as generic and reusable as possible. This approach requires slightly
424
         // in the interests of making them as generic and reusable as possible. This approach requires slightly
408
         // more multiplies, but the same number of kernel evaluations, which is the expensive part.
425
         // more multiplies, but the same number of kernel evaluations, which is the expensive part.
409
         // Inner product and scale
426
         // Inner product and scale
410
-
427
+        int ilagr = nlag-1;
411
         for (int ilag=0; ilag<nlag; ++ilag) {
428
         for (int ilag=0; ilag<nlag; ++ilag) {
412
             for (int ir2=0; ir2<nrel; ++ir2) {
429
             for (int ir2=0; ir2<nrel; ++ir2) {
413
                 //Zans(ilag, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
430
                 //Zans(ilag, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
414
-                Zans(ilag, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/Arg(ilag);
431
+                //Zans(ilag, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder()+1))/Arg(ilag);
432
+                // Segment
433
+                //std::cout << Zwork.col(ir2).segment(ilag,201).transpose() << std::endl;;
434
+                //WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder()+1).segment(ilag,201);// / Arg(ilag);
435
+                Zans(ilagr, ir2) = Zwork.col(ir2).segment(ilag,201).dot( WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder()+1) ) / Arg(ilagr);
415
             }
436
             }
437
+            ilagr -= 1;
416
         }
438
         }
439
+        //std::cout << "Zans" << Zans << std::endl;
440
+        //exit(EXIT_SUCCESS);
417
 
441
 
418
         // make sure vectors are empty
442
         // make sure vectors are empty
419
         splineVecReal.clear();
443
         splineVecReal.clear();
422
         // Now do cubic spline
446
         // Now do cubic spline
423
         // TODO Check that knots are set in right order, Eigen has reverse()
447
         // TODO Check that knots are set in right order, Eigen has reverse()
424
         //std::cout << "Arg\n" << Arg << std::endl;
448
         //std::cout << "Arg\n" << Arg << std::endl;
425
-        //std::cout << "Zans\n" << Zans.col(0) << std::endl;
449
+        //std::cout << "Arg.reverse()\n" << Arg.reverse() << std::endl;
450
+        //VectorXr Argr = Arg.reverse();
451
+        //std::cout << "Zans\n" << Zans << std::endl;
452
+        //exit(EXIT_SUCCESS);
426
         for (int ii=0; ii<Zans.cols(); ++ii) {
453
         for (int ii=0; ii<Zans.cols(); ++ii) {
427
-            auto Spline = CubicSplineInterpolator::NewSP();
428
-            Spline->SetKnots( Arg, Zans.col(ii).real() );
429
-            splineVecReal.push_back(Spline);
454
+            auto SplineR = CubicSplineInterpolator::NewSP();
455
+            SplineR->SetKnots( Arg, Zans.col(ii).real() );
456
+            splineVecReal.push_back(SplineR);
430
 
457
 
431
             auto SplineI = CubicSplineInterpolator::NewSP();
458
             auto SplineI = CubicSplineInterpolator::NewSP();
432
             SplineI->SetKnots( Arg, Zans.col(ii).imag() );
459
             SplineI->SetKnots( Arg, Zans.col(ii).imag() );
433
             splineVecImag.push_back(SplineI);
460
             splineVecImag.push_back(SplineI);
434
         }
461
         }
435
-
436
         return ;
462
         return ;
437
     }		// -----  end of method FHTKey201::ComputeLaggedRelated  -----
463
     }		// -----  end of method FHTKey201::ComputeLaggedRelated  -----
438
 
464
 

Loading…
Cancel
Save