#include "nuosc.h"
#include "earth.h"
#include "NuTransODE.h"
#include <iostream>

int main (int argc, char *argv[])
{
    using namespace blitz;

    if (argc != 9) {
        cerr << "usage: " << argv[0] 
             << " baseline"
             << " dm2_sol dm2_atm theta12 theta23 theta13 cp_phase"
             << " use_matter_effects?\n"
             << "baseline in km, dm2 in eV^2, angles in degrees\n"
             << endl;
        exit (1);
    }
        
    double baseline = atof(argv[1])*1e5;
    double dm2_sol = atof(argv[2]);
    double dm2_atm = atof(argv[3]);
    double theta12 = atof(argv[4])*M_PI/180.0;
    double theta23 = atof(argv[5])*M_PI/180.0;
    double theta13 = atof(argv[6])*M_PI/180.0;
    double cpphase = atof(argv[7])*M_PI/180.0;
    bool use_matter_effects = (argv[8][0] != '0');

    ComplexMatrix u(Range(1,3),Range(1,3));

    u = mixing_matrix(theta12,theta23,theta13,cpphase);

    cerr << "Using U=" << u; 
    cerr << "With" << (use_matter_effects?"":" no") << " matter effects\n";
    
#if 0
    NuTransODE ode(u,baseline,dm2_sol,dm2_atm,1.0e9);
    ComplexVector nu(Range(1,3));
    nu = 0;
    nu(2) = 1.0;                // nu_mu
    double step_size = 0.001*baseline;
    for (double x = 0; x<=baseline; x += step_size) {
        nu = ode.Step(x,step_size, nu);
        cout << x << " " 
             << abs(nu(1)*complex<double>(real(nu(1)),-imag(nu(1)))) << " "
             << abs(nu(2)*complex<double>(real(nu(2)),-imag(nu(2)))) << " "
             << abs(nu(3)*complex<double>(real(nu(3)),-imag(nu(3)))) << " "
             << endl;
    }

#else

    const double precision = 0.005;

    double x0[20], xf[20];
    int nregions = earth_get_slant_distances(x0,xf,baseline);
    double step_size[20];

    for (int n=0; n < nregions; ++n) {
        step_size[n] = fabs(xf[n]-x0[n])*precision;

        cerr << "    (" << x0[n] << " <--> " << xf[n] << ")\n";
        cerr << "    " 
             << earth_density(x0[n],baseline) << " -- " 
                 << earth_density(xf[n],baseline) << endl;


    }


    for (double logenergy = 8; logenergy <= 11; logenergy += 0.001) {
        double energy = pow(10.0,logenergy);
        NuTransODE ode(u,baseline,dm2_sol,dm2_atm,energy,use_matter_effects);
//        ode.SetStepper(euler_stepper);
        ComplexVector nu(Range(1,3));
        ComplexVector start(Range(1,3));
        start = 0;
        start(1) = 1;           // nu_e beam
//        nu = ode.Solve(0,baseline,baseline*precision,start);
        nu = ode.Solve(x0,xf,step_size,nregions,start);
        cout << energy << " " 
             << abs(nu(1)*complex_conjugate(nu(1))) << " "
             << abs(nu(2)*complex_conjugate(nu(2))) << " "
             << abs(nu(3)*complex_conjugate(nu(3))) << " "
//             << nu(1) << " " 
//             << nu(2) << " " 
//             << nu(3) << " "
             << vector_magnitude(nu) << endl;

    }
#endif

    return 0;
} // end of main()