// Independent probes of the unchanged upstream puncture/source provider.
#include "EffectiveSource-acceleration.h"
#include <gsl/gsl_sf_legendre.h>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>

int main() {
  std::cout << std::setprecision(17);
  for (double r0 : {6., 10.}) {
    double E=(1.-2./r0)/std::sqrt(1.-3./r0);
    double L=r0/std::sqrt(r0-3.);
    // Average radial derivative at the worldline must sum to B_r. The
    // derivative jump must supply the A_r^+-A_r^- term.
    for (int l=0;l<=16;++l) {
      double avg=0.,jump=0.;
      for (int m=l%2;m<=l;m+=2) {
        EffectiveSource source(l,m,1.);
        source.set_particle(r0,0.,0.,E,L,0.,0.,0.,0.,0.,0.);
        source.set_time_window(1.,0.,0.);
        auto field=source.Phi(r0)/r0;
        auto plus=source.dPhi_dr(r0,1)/r0-field/r0;
        auto minus=source.dPhi_dr(r0,-1)/r0-field/r0;
        double weight=(m==0?1.:2.)*gsl_sf_legendre_sphPlm(l,m,0.);
        avg+=weight*std::real((plus+minus)/2.);
        jump+=weight*std::real(plus-minus);
      }
      std::cout << "AB," << r0 << ',' << l << ',' << avg << ',' << jump << '\n';
    }
    // Approach the particle from both sides for representative modes. This
    // diagnoses regularity of each harmonic source, not of the unsummed 3D field.
    for (int l : {0,2,8,16}) {
      int m=l;
      EffectiveSource source(l,m,1.);
      source.set_particle(r0,0.,0.,E,L,0.,0.,0.,0.,0.,0.);
      source.set_time_window(1.,0.,0.);
      for(double h : {0.,1.e-2,1.e-3,1.e-4,1.e-5,1.e-6}) {
        double radii[2]={r0-h,r0+h},win[2]={1.,1.},zero[2]={0.,0.};
        double sr[2],si[2];
        source(2,radii,win,zero,zero,sr,si);
        std::cout << "source," << r0 << ',' << l << ',' << h << ','
                  << sr[0] << ',' << sr[1] << ',' << si[0] << ',' << si[1] << '\n';
      }
      for (double offset : {-1.,-.2,.2,1.}) {
        double r=r0+offset, h=1.e-4;
        int sign=offset>0?1:-1;
        double W=1.,zero=0.,sr,si;
        source(1,&r,&W,&zero,&zero,&sr,&si);
        auto second=(-source.dPhi_dr(r+2*h,sign)+8.*source.dPhi_dr(r+h,sign)
                     -8.*source.dPhi_dr(r-h,sign)+source.dPhi_dr(r-2*h,sign))/(12.*h);
        double f=1.-2./r, V=f*(l*(l+1.)/(r*r)+2./(r*r*r));
        double omega=m/std::pow(r0,1.5);
        auto independent=(omega*omega-V)*source.Phi(r)
                         +f*2./(r*r)*source.dPhi_dr(r,sign)+f*f*second;
        std::cout << "operator," << r0 << ',' << l << ',' << offset << ','
                  << std::abs(independent-std::complex<double>(sr,si)) << ','
                  << std::abs(std::complex<double>(sr,si)) << '\n';
      }
    }
  }
}