tests_data_and_calc.cpp

//
//  tests_data_and_calc.cpp
//  
//
//  Breakout by Heather Ratcliffe on 27/02/2017.
//
//

#ifdef RUN_TESTS_AND_EXIT

#include <stdio.h>
#include <math.h>
#include <cmath>
#include <mpi.h>
#include <functional>
#include "tests_data_and_calc.h"
#include "reader.h"
#include "plasma.h"
#include "controller.h"
#include "data_array.h"
#include "spectrum.h"
#include "d_coeff.h"
#include "non_thermal.h"

#include <math.h>
#include <boost/math/special_functions.hpp>
#include "valgrind.h"

test_entity_plasma::test_entity_plasma(){
/** \brief Setup tests for plasma
*
*Create plasma object from file
*/
  name = "plasma";
  //Have to set-up deck constants before creating plasma
  get_deck_constants(tests_src_dir + "test");
  share_consts();

  plas = new plasma(tests_src_dir + "test");

}
test_entity_plasma::~test_entity_plasma(){
/** \brief Teardown plasma tests*/
  delete plas;
}

int test_entity_plasma::run(){
/** \brief Test resonant frequencies and refractive indices
*
*Checks the resonant frequencies obey the equations used to derive them. Checks the dispersion roots for Whistlers match those found using high-density approx. Checks plasma O and X mode dispersion too. Note, first call for issues with these tests is to check returned mu.err on failing tests
@return Error code
*/

  int err = TEST_PASSED;
  get_deck_constants(tests_src_dir + "test");
  share_consts();
  if(!plas->is_good()){
    err |= TEST_ASSERT_FAIL;
    return err;
  }

  err |= analytic_dispersion();
  err |= high_density();
  err |= resonant_freq();
  err |= resonant_freq_poly();
  err |= other_modes();
  err |= phi_dom();

  return err;
}

int test_entity_plasma::analytic_dispersion(){
/** \brief Check analytic dispersion relations
*
*Checks the analytic relations, both ways and including derivatives
@return Error code
*/

  int err = TEST_PASSED;

  calc_type k, om, om_new, d_om, d_k;
  calc_type om_ce = plas->get_omega_ref("ce"), om_pe = plas->get_omega_ref("pe");
  size_t n_tests = 5;
  //First whistlers
  //At very large k we should get om_ce and at 0, we get 0
  k = 100;
  om = plas->get_dispersion(k, WAVE_WHISTLER);
  if(std::abs(om /om_ce -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  k = 0;
  om = plas->get_dispersion(k, WAVE_WHISTLER);
  if(std::abs(om) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  //Symmetry
  calc_type om_2;
  k = 0.0001;
  om = plas->get_dispersion(k, WAVE_WHISTLER);
  k = -k;
  om_2 = plas->get_dispersion(k, WAVE_WHISTLER);
  
  if(std::abs(om - om_2) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  
  //If we pick a few e.g.s and do them both ways we should get same result back. Start with omega for simplicity
  for(size_t i=1; i<= n_tests; i++){
    om = (float)i * om_ce/(float) (n_tests +1);
    k = plas->get_dispersion(om, WAVE_WHISTLER, true);
    om_new = plas->get_dispersion(k, WAVE_WHISTLER);
    if(std::abs(om /om_new -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  }
  
  //Derivs should be 0 at 0 for all cases and all modes
  k=0;
  //Cheat and loop through the wave indices
  for(int wave_id = WAVE_WHISTLER; wave_id <= WAVE_O; wave_id++){
    d_om = plas->get_dispersion(k, wave_id, false, true);
    if(std::abs(d_om) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  }
  //For whistler at large omega deriv is 0
  k = 100;
  d_om = plas->get_dispersion(k, WAVE_WHISTLER, false, true);
  if(std::abs(d_om) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  
  //Finally derivs at our samples should be inverses
  for(size_t i=1; i<= n_tests; i++){
    om = (float)i * om_ce/(float) (n_tests +1);
    k = plas->get_dispersion(om, WAVE_WHISTLER, true);
    
    d_om = plas->get_dispersion(k, WAVE_WHISTLER, false, true);
    d_k = plas->get_dispersion(om, WAVE_WHISTLER, true, true);
    if(std::abs(d_om*d_k - 1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  }
  //Quick check that omega is propto cos theta for fixed k, as is the derivative
  my_type theta, ref_om, ref_d_om;
  for(size_t i = 0; i< n_tests; i++){
    //Some rough range of k's
    k = 0.001;
    theta = (float) i / (float) n_tests * ANG_MAX;
    om = plas->get_dispersion(k, WAVE_WHISTLER, false, false, theta);
    d_om = plas->get_dispersion(k, WAVE_WHISTLER, false, true, theta);
    if(i == 0){
      ref_om = om;
      ref_d_om = d_om;
    }
    if(std::abs(om -ref_om * cos(theta))/om > GEN_PRECISION || std::abs(d_om - ref_d_om * cos(theta))/d_om > GEN_PRECISION){
      err |= TEST_WRONG_RESULT;
      test_bed -> report_info("Wrong angle behaviour of Whistler dispersion", 2);
    }
  }
  
  //Finally we check how it handles out of range
  try{
    k = plas->get_dispersion(om_ce*1.5, WAVE_WHISTLER, true);
  }catch(const std::exception& e){
    std::string message = e.what();
    test_bed->report_info("Exception in analytic dispersion, message " +message, 1);
    err |= TEST_ASSERT_FAIL;
  }

  //Now EM and plasma
  //At 0, we get om_pe for P and O, and X-mode cuts for X
  k = 0;
  om = plas->get_dispersion(k, WAVE_O);
  if(std::abs(om - om_pe) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  om = plas->get_dispersion(k, WAVE_PLASMA);
  if(std::abs(om - om_pe) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  om = plas->get_dispersion(k, WAVE_X_LOW);
  my_type X_cut = 0.5 * (- om_ce +  std::sqrt(om_ce*om_ce + 4.0*om_pe*om_pe));
  if(std::abs(om - X_cut) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  om = plas->get_dispersion(k, WAVE_X_UP);
  X_cut = 0.5 * (om_ce +  std::sqrt(om_ce*om_ce + 4.0*om_pe*om_pe));
  if(std::abs(om - X_cut) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  //And back again
  X_cut = 0.5 * (- om_ce +  std::sqrt(om_ce*om_ce + 4.0*om_pe*om_pe))*(1.0 + GEN_PRECISION/5.0);//Tiny bump to ensure we're in solvable range
  k = plas->get_dispersion(X_cut, WAVE_X_LOW, true);
  if(std::abs(k - 0.0) > GEN_PRECISION) err |= TEST_WRONG_RESULT;
  X_cut = 0.5 * ( om_ce +  std::sqrt(om_ce*om_ce + 4.0*om_pe*om_pe))*(1.0 + GEN_PRECISION/5.0);
  k = plas->get_dispersion(X_cut, WAVE_X_UP, true);
  if(std::abs(k - 0.0) > GEN_PRECISION) err |= TEST_WRONG_RESULT;
  
  //If we pick a few e.g.s and do them both ways we should get same result back. Go from say om_pe to 3om_pe
  for(size_t i=1; i<= n_tests; i++){
    om = plas->get_omega_ref("pe")*(1.0+ 2.0*(float)i/(float) (n_tests +1));
    k = plas->get_dispersion(om, WAVE_O, true);
    om_new = plas->get_dispersion(k, WAVE_O);
    if(std::abs(om /om_new -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
    k = plas->get_dispersion(om, WAVE_PLASMA, true);
    om_new = plas->get_dispersion(k, WAVE_PLASMA);
    if(std::abs(om /om_new -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
/*    k = plas->get_dispersion(om, WAVE_X_LOW, true);
    om_new = plas->get_dispersion(k, WAVE_X_LOW);
    if(std::abs(om /om_new -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;*/
    k = plas->get_dispersion(om, WAVE_X_UP, true);
    om_new = plas->get_dispersion(k, WAVE_X_UP);
    if(om > X_cut){
      if(std::abs(om /om_new -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
    }
  }

  //Plasma wave goes to sqrt(3) v_t and O, X to c
  k = 100;
  for(int wave_id = WAVE_O; wave_id <= WAVE_O; wave_id++){
    d_om = plas->get_dispersion(k, wave_id, false, true);
    if(std::abs(d_om - v0)/v0 > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  }
  d_om = plas->get_dispersion(k, WAVE_PLASMA, false, true);
  if(std::abs(d_om - std::sqrt(3)*0.01*v0)/d_om > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  //Finally derivs at our samples should be inverses
  for(size_t i=1; i<= n_tests; i++){
    om = om_pe + (float)i *om_pe/(float) (n_tests +1);
    for(int j = WAVE_PLASMA; j<= WAVE_O; j++){
      k = plas->get_dispersion(om, j, true);
      d_om = plas->get_dispersion(k, j, false, true);
      d_k = plas->get_dispersion(om, j, true, true);
      if(std::abs(d_om*d_k - 1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
    }
  }

  //Finally we check how it handles out of range
  try{
    for(int j = WAVE_PLASMA; j < WAVE_X_LOW; j++){
      k = plas->get_dispersion(plas->get_omega_ref("pe")*0.3, j, true);
    }
  }catch(const std::exception& e){
    std::string message = e.what();
    test_bed->report_info("Exception in analytic dispersion, message " +message, 1);
    err |= TEST_ASSERT_FAIL;
  }

  if(err == TEST_PASSED) test_bed->report_info("Analytic dispersion OK", 1);
  else  test_bed->report_info("Error in analytic dispersion", 1);

  return err;
}

int test_entity_plasma::resonant_freq(){
/** \brief Check resonant frequency solver
*
*Checks the returned resonant frequency obeys equations used to derive it by solving both for mu. 
@return Error code
*/

  int err=TEST_PASSED;
  int n_tests = 10;
  std::vector<calc_type> results;
  mu_dmudom my_mu;
  int err_count = 0;
  calc_type x, v_par, n, om_ce_local, om_pe_local, om_ce_signed;
  om_ce_local = plas->get_omega_ref("ce");
  om_pe_local = plas->get_omega_ref("pe");
  om_ce_signed = - std::abs(om_ce_local);
  calc_type cos_theta, mu_tmp1, mu_tmp2;
  calc_type gamma;

  test_bed->report_info("Testing resonant frequency solver", 1);

  //Check the n=0 and v=0 degenerate cases
  //Zero angle, zero velocity, any n, expect empty result
  results = plas->get_resonant_omega(0.0, 0.0, 1.0, -1);
  if(results.size() != 1 && (results[0]/om_ce_local - 1.0 > PRECISION)){
    test_bed->report_info("Erroneous solution when v=0 in resonant frequency solver", 2);
    err |= TEST_WRONG_RESULT;
  }
  
  calc_type expected_result =  0.50, corresponding_v =  0.162251 * v0, omega_solution;
  gamma = gamma_rel(corresponding_v);
  results = plas->get_resonant_omega(0.0, corresponding_v, gamma, -1);
  //Insert a known result here from the IDL code, with one root at 0 angle
  //Since we're using a sample like this, only expect a few sf of equality as the IDL uses quite different process
  if(results.size() != 1 || std::abs(std::abs(results[0]/om_ce_local)-expected_result) > LOW_PRECISION){
    test_bed->report_info("Erroneous solution for example case in resonant frequency solver", 2);
    err |= TEST_WRONG_RESULT;
  }
  //Check over the range of other cases
  //Loop over particle velocity, assuming propagation at alpha = pi/8 say
  calc_type theta, test_angle = pi/8.0;
  for(int ii=0; ii<n_tests; ii++){
    v_par = (0.01 + 0.5*(float)ii/ (float)(n_tests+1))* v0;
    //Loop over angles
    for(int j=0; j< n_tests; j++){
      x = 4.0 * (float) j / (float)(n_tests+1);
      theta = std::atan(x);
      cos_theta = std::cos(theta);
      //Loop over n
      for(int k=0; k< n_tests; k++){
        n = -n_tests/2 + k*n_tests/2;
        
        gamma = gamma_rel(v_par/cos(test_angle));
        
        results = plas->get_resonant_omega(theta, v_par, gamma, n);
        /**Now check each element of the resonant frequency solution set satisfies Stix 2.45 and the resonance condition together*/
        for(size_t i=0; i<results.size(); ++i){
          omega_solution = results[i];

          //Solve Res condition for mu^2 = (kc/om)^2
          mu_tmp1 = std::pow(v0 * (gamma * omega_solution - n*om_ce_signed)/(gamma * omega_solution * v_par *cos_theta), 2);
          mu_tmp2 = (1.0 - (std::pow(om_pe_local,2)/(omega_solution *(omega_solution + om_ce_signed*cos_theta))));
          //Since we're solving a cubic, expect roundoff error to be exacerbated, so we bump the threshold to roughly (5e-16)^(1/3) (cube root epsilon) and I've tightened by a factor 10 also since I'm doing fractional
          if(std::abs((mu_tmp1 - mu_tmp2)/mu_tmp1) > 2.5*NUM_PRECISION){
            err|=TEST_WRONG_RESULT;
            test_bed->report_info("Refractive index mismatch of "+mk_str((std::abs((mu_tmp1-mu_tmp2))/mu_tmp1)*100) +'%', 2);
          }

          //Also check against the full-mu solution
          my_mu = plas->get_high_dens_mu(omega_solution, theta);
          if(my_mu.err){
            err|=TEST_WRONG_RESULT;
            test_bed->report_info("No full mu solution for resonant frequency", 2);
            err_count++;
          }else{
            //This should match the high-density mu solution almost exactly, once we remove the leading 1.0
            if(std::abs(my_mu.mu - std::sqrt(mu_tmp2 - 1.0))/std::sqrt(mu_tmp2 - 1.0) > NUM_PRECISION){
              err|=TEST_WRONG_RESULT;
              test_bed->report_info("Mismatched full mu solution in resonant solver, error " + mk_str(std::abs(my_mu.mu-std::sqrt(mu_tmp2-1.0))/std::sqrt(mu_tmp2-1.0)*100.0, true)+'%', 2);
            }
          }
        }
      }
    }
  }
  test_bed->report_info("Mu error count: "+mk_str(err_count), 2);
  return err;

}

int test_entity_plasma::resonant_freq_poly(){
/** \brief Check resonant frequency solver
*
*Checks the returned resonant frequencies from the high-density and full polynomials have full solutions and match to reasonable degree
@return Error code
*/

  int err=TEST_PASSED;
  int n_tests = 10;
  std::vector<calc_type> results, HD_results;
  mu_dmudom my_mu;
  int err_count = 0;
  calc_type x, v_par, n, om_ce_local, om_pe_local, om_ce_signed;
  om_ce_local = plas->get_omega_ref("ce");
  om_pe_local = plas->get_omega_ref("pe");
  om_ce_signed = - std::abs(om_ce_local);
  calc_type cos_theta, mu_tmp1, mu_tmp2;
  calc_type gamma, tmp_result;

  test_bed->report_info("Testing resonant frequency polynomial solver", 1);

  //Check the n=0 and v=0 degenerate cases
  //Zero angle, zero velocity, any n, expect empty result
  results = plas->get_resonant_omega_full(0.0, 0.0, 1.0, -1);
  if(results.size() != 1 && (results[0]/om_ce_local - 1.0 > PRECISION)){
    test_bed->report_info("Erroneous solution when v=0 in resonant frequency solver", 2);
    err |= TEST_WRONG_RESULT;
  }
  
  //Check over the range of cases
  //Loop over particle velocity, assuming propagation at alpha = pi/8 say
  calc_type theta, test_angle = pi/8.0;
  for(int ii=0; ii<n_tests; ii++){
    v_par = (0.01 + 0.5*(float)ii/ (float)(n_tests+1))* v0;
    gamma = gamma_rel(v_par/cos(test_angle));
    //Loop over angles
    for(int j=0; j< n_tests; j++){
      x = 4.0 * (float) j / (float)(n_tests+1);
      theta = std::atan(x);
      cos_theta = std::cos(theta);
      //Loop over n
      for(int k = 0; k < n_tests; k++){
        n = -5 + k*10/n_tests;
        results = plas->get_resonant_omega_full(theta, v_par, gamma, n);
        HD_results = plas->get_resonant_omega(theta, v_par, gamma, n);
        //Remove any HD results outside the NR bounds
        for(size_t nn = 0; nn < HD_results.size(); nn++){
          if(HD_results[nn] > om_ce_local*NR_max_om || HD_results[nn] < om_ce_local*NR_min_om){
            HD_results.erase(HD_results.begin() + nn);
            nn--;
          }
        }
        //This check just makes sure the solver is not mis-reporting solutions, by checking they are actually solutions to the resonant poly.
        for(size_t i=0; i < results.size(); ++i){
          if(!plas->check_resonant_omega_full(theta, v_par, gamma, n, results[i]/om_ce_local, tmp_result)){
            test_bed->report_info("Non-root solution, value "+mk_str(results[i]));
            err |= TEST_WRONG_RESULT;
          }
        }

        if(results.size() != HD_results.size()){
          err |= TEST_WRONG_RESULT;
          test_bed->report_info("Mismatched solution numbers "+mk_str(results.size()) +" vs "+mk_str(HD_results.size()));
        }else{
          if(HD_results.size() > 1 && HD_results[0] > HD_results[1]){
            //Reverse cubic solutions to ensure smaller is first
            double tmp = HD_results[0];
            HD_results[0] = HD_results[1];
            HD_results[1] = tmp;
          }
          if(results.size() > 1 && results[0] > results[1]){
            //Reverse solutions to ensure smaller is first
            double tmp = results[0];
            results[0] = results[1];
            results[1] = tmp;
          }

          for(size_t i=0; i < results.size(); ++i){
            if(std::abs(results[i] - HD_results[i])/results[i] > 0.1){
              err_count++;
              err |= TEST_WRONG_RESULT;
              test_bed->report_info("HD and full solution do not match, error "+mk_str((results[i] - HD_results[i])/results[i], true) + " at "+mk_str(i), 2);
            }
          }
        }
      }
    }
  }
  test_bed->report_info("Mu error count: "+mk_str(err_count), 2);
  return err;

}

int test_entity_plasma::high_density(){
/** \brief Tests high density approximation for dispersion relations
*
*Test if the mu found by get_root and get_phi_mu_om matches the high density whistler in high dens regime
@return Error code
*/

  int err=TEST_PASSED;

  calc_type om_ce_local, om_pe_local;
  calc_type mu_tmp2;

  om_ce_local = plas->get_omega_ref("ce");
  om_pe_local = plas->get_omega_ref("pe");

  size_t n_tests = 10;
  calc_type tmp_omega = 0.0, tmp_theta = pi/(calc_type)(n_tests+1), gamma_particle = 1.0;
  mu_dmudom my_mu, my_mu_dens;
  int err_cnt = 0;
  test_bed->report_info("Testing whistler high density approx.", 1);

  for(size_t i =0; i<n_tests; i++){
    tmp_omega += std::abs(om_ce_local)/(calc_type)(n_tests + 1);
    my_mu = plas->get_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle);
    my_mu_dens = plas->get_high_dens_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle);
    /** my_mu.mu should roughly equal Stix 2.45*/
    mu_tmp2 = sqrt(1.0 - (std::pow(om_pe_local,2)/(tmp_omega*(tmp_omega + om_ce_local*std::cos(tmp_theta)))));
    if(std::abs(my_mu.mu-mu_tmp2)/my_mu.mu > LOW_PRECISION){
      err_cnt++;
      test_bed->report_info("Mismatch in high density approx or dispersion solver at omega="+mk_str(tmp_omega/std::abs(om_ce_local), true)+" om_ce and theta= "+mk_str(tmp_theta/pi)+" pi", 1);
      test_bed->report_info("Mu "+mk_str(my_mu.mu)+" difference "+mk_str(my_mu.mu - mu_tmp2)+" relative error "+mk_str((my_mu.mu-mu_tmp2)/my_mu.mu, true), 2);
    }
    /** my_mu_dens should EXACTLY equal Stix 2.45 without the 1.0 term*/
    mu_tmp2 = sqrt( - (std::pow(om_pe_local,2)/(tmp_omega*(tmp_omega + om_ce_local*std::cos(tmp_theta)))));

    if(std::abs(my_mu_dens.mu-mu_tmp2)/my_mu_dens.mu > NUM_PRECISION){
      err_cnt++;
      test_bed->report_info("  Mismatch in high density approx or dispersion solver at omega="+mk_str(tmp_omega/std::abs(om_ce_local), true)+" om_ce and theta= "+mk_str(tmp_theta/pi)+" pi", 1);

      test_bed->report_info("  Mu "+mk_str(my_mu_dens.mu)+" difference "+mk_str(my_mu_dens.mu - mu_tmp2)+" relative error "+mk_str((my_mu_dens.mu-mu_tmp2)/my_mu_dens.mu), 2);
    }
  }
  
  tmp_omega = 0.6*std::abs(om_ce_local);
  for(size_t i =0; i<n_tests; i++){
    tmp_theta += pi/(calc_type)(n_tests+1);
    my_mu = plas->get_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle);
    /* my_mu.mu should roughly equal Stix 2.45*/
    mu_tmp2 = sqrt(1.0 - (std::pow(om_pe_local,2)/(tmp_omega*(tmp_omega + om_ce_local*std::cos(tmp_theta)))));
    if(std::abs(my_mu.mu-mu_tmp2)/my_mu.mu > LOW_PRECISION){
      err_cnt++;
      test_bed->report_info("    Mismatch in high density approx or dispersion solver at omega="+mk_str(tmp_omega/std::abs(om_ce_local), true)+" om_ce and theta= "+mk_str(tmp_theta/pi)+" pi", 1);

      test_bed->report_info("    Mu "+mk_str(my_mu.mu)+" difference "+mk_str(my_mu.mu - mu_tmp2)+" relative error "+mk_str((my_mu.mu-mu_tmp2)/my_mu.mu), 2);
    }
  }
  if(err_cnt> 0){
    test_bed->report_info("Total "+mk_str(err_cnt)+" out of "+mk_str(2*(int)n_tests)+" issues in high density approx or dispersion solver at precision: "+mk_str(LOW_PRECISION), 1);
    //err|=TEST_WRONG_RESULT;
    //Make these a warning not an error because we expect them sometimes
  }

  return err;
}

int test_entity_plasma::other_modes(){
/** \brief Test dispersion solver further
*
*Tests dispersion solver for other wave modes, O and X
@return Error code
*/

  int err = TEST_PASSED;

  calc_type om_ce_local, om_pe_local;
  calc_type mu_tmp2, gamma_particle = 1.0;

  om_ce_local = plas->get_omega_ref("ce");
  om_pe_local = plas->get_omega_ref("pe");


  test_bed->report_info("Testing dispersion solver for plasma O mode", 1);
  size_t n_tests = 10;
  mu_dmudom my_mu;
  /*Try plasma wave modes in solvers, perpendicular propagation*/
  calc_type tmp_omega = om_pe_local;
  calc_type tmp_theta = pi/2.0;
  for(size_t i =0; i<n_tests; i++){
    tmp_omega += std::abs(om_pe_local)/(calc_type)(n_tests + 1);
    my_mu = plas->get_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle);

    mu_tmp2 = std::sqrt(std::pow(tmp_omega, 2) - std::pow(om_pe_local, 2))/tmp_omega;
    
    if(std::abs(my_mu.mu-mu_tmp2)/my_mu.mu > LOW_PRECISION){
      test_bed->report_info("Error in approx or dispersion solver for plasma wave at "+mk_str(tmp_omega/std::abs(om_pe_local))+" "+mk_str(tmp_theta), 1);
      test_bed->report_info("Mu "+mk_str(my_mu.mu)+" difference "+mk_str(my_mu.mu - mu_tmp2)+" relative error "+mk_str((my_mu.mu-mu_tmp2)/my_mu.mu), 2);
    }
  }
  /*Try left hand X mode too*/
  test_bed->report_info("Testing dispersion solver for plasma X mode", 1);

  calc_type omega_UH = std::sqrt(om_pe_local*om_pe_local + om_ce_local*om_ce_local);
  for(size_t i = 0; i < n_tests; i++){
    tmp_omega += std::abs(om_pe_local)/(calc_type)(n_tests + 1);
    my_mu = plas->get_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle, false);
    mu_tmp2 = std::sqrt(1.0 - std::pow(om_pe_local/tmp_omega, 2)*(std::pow(tmp_omega, 2) - std::pow(om_pe_local, 2))/(std::pow(tmp_omega, 2) - std::pow(omega_UH, 2)));

    if(std::abs(my_mu.mu-mu_tmp2)/my_mu.mu > LOW_PRECISION){
      test_bed->report_info("Error in approx or dispersion solver for plasma wave at "+mk_str(tmp_omega/std::abs(om_pe_local))+" "+mk_str(tmp_theta), 1);
      test_bed->report_info("Mu "+mk_str(my_mu.mu)+" difference "+mk_str(my_mu.mu - mu_tmp2)+" relative error "+mk_str((my_mu.mu-mu_tmp2)/my_mu.mu), 2);
    }
  }
  return err;
}

int test_entity_plasma::phi_dom(){
/** \brief Test other plasma returns
*
* Checks the values of mu.dom, mu.dmudtheta and phi against special cases.
@return Error code
*/

  int err = TEST_PASSED;
  size_t n_tests = 10;

  calc_type mu_tmp1, om_ce_local, tmp_phi;
  om_ce_local = plas->get_omega_ref("ce");

  calc_type d_omega = std::abs(om_ce_local)/1e8;
  calc_type d_theta = pi/(calc_type)(n_tests)/1e7;
  //Derivative step size.

  mu_dmudom my_mu, my_mu_p;

  calc_type tmp_omega = 0.0, tmp_theta = pi/(calc_type)(n_tests), gamma_particle = 1.0;

  test_bed->report_info("Testing phi", 1);
  //For phi all we can really do is check that it's positive (as it's a square) and check it's not NaN, especially for n=0, and also of "order 1"
  for( int n = -4; n < 4 ; n++){
    my_mu = plas->get_high_dens_phi_mu_om(std::abs(om_ce_local)* 0.5, tmp_theta, pi/10.0, 0, gamma_particle);
    if(my_mu.phi != my_mu.phi || my_mu.phi < 0){
      err |= TEST_ASSERT_FAIL;
      test_bed->report_info("Phi is invalid in get_high_dens_phi_mu_om for n="+mk_str(n), 2);
    }
    tmp_phi = my_mu.phi;
    my_mu = plas->get_phi_mu_om(std::abs(om_ce_local)* 0.5, tmp_theta, pi/10.0, 0, gamma_particle);
    if(my_mu.phi != my_mu.phi || my_mu.phi < 0){
      err |= TEST_ASSERT_FAIL;
      test_bed->report_info("Phi is invalid in get_phi_mu_om for n="+mk_str(n), 2);
    }
    //Allow 15%
    if((my_mu.phi- tmp_phi)/my_mu.phi > 0.15){
      err |= TEST_WRONG_RESULT;
      test_bed->report_info("Inconsistent phi ("+mk_str(my_mu.phi)+ ") with high dens phi ("+mk_str(tmp_phi)+") for n="+mk_str(n), 2);
    }
    if(tmp_phi < 0.1 || tmp_phi > 10.0){
      err |= TEST_WRONG_RESULT;
      test_bed->report_info("Phi is not order 1 for n="+mk_str(n), 2);
    }
  }

  test_bed->report_info("Testing dmu/domega", 1);

  for(size_t i =0; i<n_tests; i++){
    tmp_omega += std::abs(om_ce_local)/(calc_type)(n_tests + 1);
    my_mu = plas->get_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle);
    my_mu_p = plas->get_phi_mu_om(tmp_omega+d_omega, tmp_theta, 0.0, 0, gamma_particle);

    /** Approx numerical derivative*/
    mu_tmp1 = (my_mu.mu - my_mu_p.mu)/d_omega;
    if(std::abs(std::abs(mu_tmp1/my_mu.dmudom) - 1.0) > NUM_PRECISION){
      err|=TEST_WRONG_RESULT;
      test_bed->report_info("Wrong derivative in get_phi_mu_om", 2);
    }
  }
  
  test_bed->report_info("Testing dmu/dtheta", 1);
  tmp_theta = 0.0;
  tmp_omega = 0.7*std::abs(om_ce_local);
  for(size_t i =0; i<n_tests; i++){
    tmp_theta += pi/(calc_type)(n_tests+1);
    my_mu = plas->get_phi_mu_om(tmp_omega, tmp_theta, 0.0, 0, gamma_particle);
    my_mu_p = plas->get_phi_mu_om(tmp_omega, tmp_theta+d_theta, 0.0, 0, gamma_particle);

    /** Approx numerical derivative. Manually fix signs*/
    mu_tmp1 = -(my_mu.mu - my_mu_p.mu)/d_theta;
    if(std::abs(std::abs(mu_tmp1 /my_mu.dmudtheta) - 1.0) > NUM_PRECISION){
      err|=TEST_WRONG_RESULT;
      test_bed->report_info("Wrong derivative in get_phi_mu_om at omega = "+mk_str(tmp_omega/std::abs(om_ce_local), true) +" and phi = "+mk_str(tmp_theta/pi, true)+" pi", 2);

    }
  }

  return err;

}

/*
For getting the resonant frequency to consider errors will matter, but the noise in spectra can be expected to be a similar sort of error
*
The reason for using the better dispersion solver is a) to avoid any numerical derivatives in general and b) because some of the calculation depends on derivs of the refractive index

*So basically I am using the easy approx where a) it’s a genuine nightmare to do better (10-14th order polynomial) and b) I expect other errors to be similar or more important and c) I really hope it’s linear or polynomial in error*/
//----------------------------------------------------------------

test_entity_spectrum::test_entity_spectrum(){
/** \brief Setup spectrum tests*/
  name = "spectrum checks";
  file_prefix = tests_src_dir;
}
test_entity_spectrum::~test_entity_spectrum(){
/**\brief Teardown spectrum tests*/
  if(test_contr) delete test_contr;
}

int test_entity_spectrum::run(){
/** \brief Test spectrum extraction
*
* This tests the dispersion relation approximations are OK (plain and vg). Check test spectrum makes sense. Test extraction of a spectrum from data. Note data does not come from files, but from a test file which is already written as a data array using generate_fftd.pro which makes FFT_data.dat with the FFTd data and spectrum.dat with a derived spectrum to check against.
@return Error code
*/

  int err = TEST_PASSED;

  get_deck_constants(file_prefix);
  share_consts();

  err|= setup();
  if(!test_bed->check_for_abort(err)) err|= basic_tests1();
  if(!test_bed->check_for_abort(err)) err|= basic_tests2();
  if(!test_bed->check_for_abort(err)) err|= technical_tests();
  if(!test_bed->check_for_abort(err)) err|= albertGs_tests();
  
  return err;

}

int test_entity_spectrum::setup(){
/** \brief Setup to test spectrum
*
* Read test FFT data from file, and create controller. Note strictly this is the test of data array constructor taking a filename too.
@return Error code
*/

  int err = TEST_PASSED;

  test_dat_fft = data_array(file_prefix + "FFT_data.dat");
  if(!test_dat_fft.is_good()){
    err |= TEST_ASSERT_FAIL;
    err |= TEST_FATAL_ERR;
  }
  test_contr = new controller(file_prefix);

  return err;
}

int test_entity_spectrum::basic_tests1(){
/** \brief Basic tests of spectrum
*
*Test test_spectrum and angle generation
@return Error code
*/
  int err = TEST_PASSED;

  std::fstream outfile;
  size_t len=0;
  my_type total_error =0.0;
  my_type * d_angle, * angle_data;

  test_contr->add_spectrum(test_dat_fft.get_dims(0), DEFAULT_N_ANG, true);
  //Delete and re-add to test that
  test_contr->delete_current_spectrum();
  test_contr->add_spectrum(test_dat_fft.get_dims(0), DEFAULT_N_ANG, true);

  test_contr->get_current_spectrum()->make_test_spectrum(FUNCTION_DELTA);
  //Check angle distrib integrates to 1 for each case
  //NOTE we can only do this if MIN_ANG is either 0 or is - MAX_ANG. otherwise we're into erf and bunk
  bool is_symmetric=false, is_zero = false;
  if(std::abs(ANG_MIN + ANG_MAX) < PRECISION) is_symmetric = true;
  if(std::abs(ANG_MIN) < PRECISION) is_zero = true;
  
  if(is_symmetric || is_zero){
    len = test_contr->get_current_spectrum()->get_g_dims(1);
    d_angle = (my_type *) calloc(DEFAULT_N_ANG, sizeof(my_type));
    for(size_t i=0; i<DEFAULT_N_ANG-1; ++i){
      d_angle[i] = std::abs(test_contr->get_current_spectrum()->get_ang_axis_element(i) - test_contr->get_current_spectrum()->get_ang_axis_element(i+1));
    }
    angle_data = (my_type *) malloc(len*sizeof(my_type));
    for(size_t i=0; i<len; i++){
      *(angle_data + i) = test_contr->get_current_spectrum()->get_g_element(0, i);
    }
    
    total_error = integrator(angle_data, len, d_angle);
    test_contr->get_current_spectrum()->make_test_spectrum(FUNCTION_GAUSS);
    for(size_t i=0; i<len; i++) *(angle_data + i) = test_contr->get_current_spectrum()->get_g_element(0, i);

    total_error += integrator(angle_data, len, d_angle);

    test_contr->get_current_spectrum()->make_test_spectrum(FUNCTION_ISO);
    for(size_t i=0; i<len; i++) *(angle_data + i) = test_contr->get_current_spectrum()->get_g_element(0, i);
    total_error += integrator(angle_data, len, d_angle);
    
    my_type expected = is_zero ? 1.5 : 3.0;
    //Iso always integrates to 1. Gaussian and delta are always symmetric
    if(std::abs(total_error - expected)/3.0 > NUM_PRECISION){
    
      err |= TEST_WRONG_RESULT;
      test_bed->report_info("Error in angular distribution integrals, value " + mk_str(total_error, true));
    }
    if(angle_data) free(angle_data);
    if(d_angle) free(d_angle);
  }else{
    test_bed->report_info("Cannot test assymmetric spectrum");
  }
  outfile.open(tests_tmp_dir+"spect_testy.dat", std::ios::binary|std::ios::out);
  test_contr->get_current_spectrum()->write_to_file(outfile);
  outfile.close();
  if(err == TEST_PASSED) test_bed->report_info("Test spectrum OK");
  
  return err;
}
int test_entity_spectrum::basic_tests2(){
/**\brief Test spectrum extraction
*
*Compare extracted spectrum from FFT'd data file to test file
@return Error code
*/
  int err = TEST_PASSED;

  std::fstream outfile, infile;
  size_t len = 0;
  my_type total_error = 0.0;
  test_contr->get_current_spectrum()->generate_spectrum(test_dat_fft ,10, FUNCTION_GAUSS, DEFAULT_SPECTRUM_ANG_STDDEV);

  test_spect = data_array(file_prefix + "spectrum.dat");
  if(test_spect.is_good()){
    //We ignore frequencies below say 0.05 om_ce
    const my_type * ax = test_spect.get_axis(0, len);
    int min_ind = 0;
    if(ax) min_ind = where(ax+len/2, len/2, 17588.200*0.05);
    /**Hard code min freq to match the IDL file with test data generation...*/
    
    total_error = 0.0;
    for(size_t i=0; i< len/2 - min_ind; i++){
      total_error += std::abs(test_contr->get_current_spectrum()->get_B_element(i)-test_spect.get_element(i));
    }
    for(size_t i=len/2 + min_ind; i< len; i++){
      total_error += std::abs(test_contr->get_current_spectrum()->get_B_element(i)-test_spect.get_element(i));

    }
    if(total_error > LOW_PRECISION){
      err |= TEST_WRONG_RESULT;
      test_bed->report_info("Mismatch between generated spectrum and test spectrum of "+mk_str(total_error));
    }
    /* Preserve the spectrum*/
    outfile.open(tests_tmp_dir + "spect_out.dat", std::ios::binary|std::ios::out);
    test_contr->get_current_spectrum()->write_to_file(outfile);
    outfile.close();
  }else{
    err |= TEST_ASSERT_FAIL;
  }
  if(err == TEST_PASSED) test_bed->report_info("Generate spectrum OK");
  
  data_array old_B = test_contr->get_current_spectrum()->copy_out_B();
  //Now dump to file and read back in and compare
  bool err2 = test_contr->add_spectrum(tests_tmp_dir + "spect_out.dat");
  if(err2) err |= TEST_ASSERT_FAIL;
  else{
    data_array new_B = test_contr->get_current_spectrum()->copy_out_B();
    if(!old_B.is_good() || !new_B.is_good() || (old_B != new_B)){
      test_bed->report_info("Error or Mismatch in read", 0);
      err|= TEST_WRONG_RESULT;
    }
  }
  return err;

}

void test_entity_spectrum::set_vals(spectrum& spect){
/** \brief Set B and g to sample vals
*
* Sets B and g (assuming angle-is-function) to distinguishable vals
@param spect Spectrum to set values of
*/

  if(!spect.get_g_is_angle_only()) return;

  size_t om_sz = spect.get_omega_length();
  size_t ang_sz = spect.get_angle_length();
  
  for(size_t i = 0; i < om_sz; i++){
    spect.set_B_element(i, i);
  }
  for(size_t i = 0; i < ang_sz; i++){
    spect.set_g_element(i, i*2);
  }

}

int test_entity_spectrum::technical_tests(){
/** \brief Test spectrum object
*
*Checks spectrum object copy etc
@return Error code
*/
  int err = TEST_PASSED;

  test_bed->report_info("Checking technical aspects", 2);

  //Create empty spectrum
  spectrum test_spect = spectrum(10, 10, true);
  set_vals(test_spect);

  spectrum test_spect2 = test_spect;

  if(test_spect2 != test_spect){
    err |= TEST_WRONG_RESULT;
    test_bed->report_info("Copy or equality problem", 1);
  }
  try{
    std::vector<spectrum> my_vec;
    my_vec.push_back(test_spect2);
    my_vec.push_back(test_spect2);
    if(my_vec[0] != test_spect2) err|= TEST_WRONG_RESULT;
  }catch(const std::exception& e){
    //Swallow and continue if possible,
    std::string message = e.what();
    test_bed->report_info("Exception message " +message, 1);
    err |= TEST_ASSERT_FAIL;
  }
  
  if(err == TEST_PASSED) test_bed->report_info("Technical aspects OK", 1);
  return err;

}

int test_entity_spectrum::albertGs_tests(){
/** \brief Tests of the Albert G functions
*
*Tests the calculation of G_1 and G_2 in Albert \cite Albert2005 by get_G1 and get_G2. Also tests the normalisations on the way.
@return Error code
\caveat The I(omega) calc below uses a splunged version of Lyons \cite Lyons1974b A7 which is half cold, half warm plasma. Since we're only after a sanity check here, we just use restrictive angles (where cold approx performs "better") and allow 10% mismatch.
*/
  int err = TEST_PASSED;

  calc_type om_ce_local, om_pe_local, G1, G2, G1_analytic, G2_analytic, G1_tracker = 0.0;
  om_ce_local = test_contr->get_plasma().get_omega_ref("ce");
  om_pe_local = test_contr->get_plasma().get_omega_ref("pe");

  calc_type mass_ratio = 1.0/1836.2;

  size_t n_tests = 30;
  calc_type tmp_omega = 0.0, tmp_x;

  test_contr->add_spectrum(2048, DEFAULT_N_ANG, true);

  if(!test_contr->get_current_spectrum()->is_good()){
    my_error_print("Spectrum in invalid state. Aborting", mpi_info.rank);
    err |=TEST_ASSERT_FAIL;
    err |=TEST_FATAL_ERR;
    return err;
  }

  test_contr->get_current_spectrum()->make_test_spectrum(FUNCTION_GAUSS);
    
  my_type om_min, om_max, x_min, x_max, om_peak;
  om_min = 2000.0;
  om_max = 16500.0;
  //make sure this is lower than the test spectrum axis range
  x_min = 0.0;
  x_max = 0.4;
  //Because I(omega) uses COLD plasma dispersion it really doesn't work well at "large" angles, above say 30 deg. So we cut off really harshly here
  //This might be a bit delicate with respect to the degree of mismatch, but I think the problem is just because of the cold vs warm approx and is enough for a sanity check
  
  test_contr->get_current_spectrum()->truncate_om(om_min, om_max);
  test_contr->get_current_spectrum()->truncate_x(x_min, x_max);
  om_peak = test_contr->get_current_spectrum()->get_peak_omega();

  //Now we have a test spectrum. Need to know what its normalisations should be. And what the Albert functions should resolve to.
  my_type width = 0.1*om_peak;
  
  for(size_t i = 0; i < n_tests; i++){
    tmp_omega = om_min + (float) i/(float) n_tests * (om_max-om_min);
    //Cover range from small to just below om_ce...
    G1 = get_G1(test_contr->get_current_spectrum(), tmp_omega);

    //Analytic calculations for truncated Gaussians, Albert
    
    if(tmp_omega > om_min && tmp_omega < om_max){
      G1_analytic = 2.0 / std::sqrt(pi) * std::exp( - std::pow((tmp_omega - om_peak)/width, 2));
      G1_analytic /= (boost::math::erf((om_max - om_peak)/width) +boost::math::erf((om_peak - om_min)/width));
      G1_analytic /=width;
    }else{
      G1_analytic = 0.0;
    }
    G1_tracker += G1_analytic;//Keep sum to check we're not hitting zero everywhere
    //Arbitrary precision, slightly tuned to current parameters
    if( G1_analytic > 0.0 && ((G1 != 0.0 && std::abs(G1-G1_analytic)/(G1) > 0.02)|| (G1 == 0.0 && G1_analytic != 0.0))){
      err |= TEST_WRONG_RESULT;
      test_bed->report_info("G1 does not match analytic calc, relative error = "+mk_str((std::abs(G1/G1_analytic)-1.0)*100, true)+"% at "+mk_str(tmp_omega, true), mpi_info.rank);
    }
  }
  if(G1_tracker < tiny_my_type){
    err |= TEST_ASSERT_FAIL;
    test_bed->report_info("G1 is always zero", mpi_info.rank);
  }

  size_t ang_sz = test_contr->get_current_spectrum()->get_angle_length();
  calc_type I_om = 0.0;
  size_t counter = 0, none_counter = 0;
  for(size_t j = 1; j < n_tests; j++){
    tmp_omega = (float) j/(float) (n_tests) * std::abs(om_ce_local);
    I_om = calc_I_omega(tmp_omega, test_contr->get_current_spectrum(), test_contr);

    for(size_t i = 0; i < ang_sz;i++){
      counter++;
      tmp_x = test_contr->get_current_spectrum()->get_ang_axis_element(i);

      G2 = get_G2(test_contr->get_current_spectrum(), tmp_omega, tmp_x);
      G2_analytic = 0.0;
      if((tmp_omega > om_min && tmp_omega < om_max && tmp_x > x_min && tmp_x < x_max)){
        G2_analytic = std::pow((( mass_ratio / (1.0 + mass_ratio))*om_ce_local*om_ce_local/om_pe_local/om_pe_local), 1.5);
        //Same g in num and denom so don't need to normalise
        G2_analytic *= test_contr->get_current_spectrum()->get_g_element(i);
        G2_analytic /= I_om;
        /** \todo Trace this 2! I suspect it's from the 2 I can't reproduce in Lyons A7*/
        G2_analytic *= 2.0;
      }else{
        G2_analytic = 0.0;
      }
      //Both should be "non-zero" and we allow rather large mismatch. See caveat above.
      if(G2_analytic > 1e-30 && G2 > 1e-30 && std::abs(G2/G2_analytic - 1.0) > 0.1){
        none_counter++;
        test_bed->report_info("G2 does not match analytic calc, relative error = "+mk_str((std::abs(G2/G2_analytic)-1.0)*100, true)+"% at omega="+mk_str(tmp_omega/om_ce_local, true)+" and x="+mk_str(tmp_x, true), mpi_info.rank);
        err |= TEST_WRONG_RESULT;
      }
    }
  }
  if(none_counter > 0) my_error_print("Tested "+mk_str(counter)+" and got "+mk_str(none_counter)+" errors", mpi_info.rank);
  std::fstream outfile;
  outfile.open(tests_tmp_dir + "spect_truncated.dat", std::ios::out|std::ios::binary);
  test_contr->get_current_spectrum()->write_to_file(outfile);
  outfile.close();


  return err;
}

my_type calc_I_omega(my_type omega, spectrum * my_spect, controller * my_contr){
/** \brief Calculate the function I(omega)
*
* Calculates I(omega) as in Lyons \cite Lyons1974B directly for a Gaussian g. Note I am mixing cold and warm plasma theory, so this is only a broad sanity check.
@param omega Frequency to calculate at
@param my_spect Spectrum object to use
@param my_contr Controller providing plasma
@return Value of I(omega)
*/
  my_type Psi, Psi2, theta, x, dx, g_x, I_contrib, I_contrib2, I_om = 0.0, om_sq_p_e;
  size_t x_sz = my_spect->get_angle_length();

  my_type om_ce_local = my_contr->get_plasma().get_omega_ref("ce");
  my_type om_pe_local = my_contr->get_plasma().get_omega_ref("pe");
  calc_type M = 1.0/1836.2;//m_e/m_p
  my_type om_cp = om_ce_local*M;
  om_sq_p_e = omega*omega/om_ce_local/om_cp;
  mu_dmudom my_mu;
  plasma my_plas = my_contr->get_plasma();
  I_om = 0.0;
  for(size_t i = 1; i < x_sz; i++){
    x = TAN_MAX * (float)i/ (float)x_sz;
    dx = TAN_MAX/(float)x_sz;
    theta = std::atan(x);

    my_mu = my_plas.get_mu(omega, theta);
    
    if(!my_mu.err && std::abs(omega) < std::abs(om_ce_local*cos(theta))){
      //Mu has no solutions where omega exceeds Om_ce*cos(theta), or when there's an err
      //NB sign selected for Whistler branch
      //Psi from cold plasma theory
      Psi2 = 1.0 - om_sq_p_e - std::pow(std::sin(theta), 2)/2.0 + std::sqrt(std::pow(std::sin(theta), 4)/4.0 + std::pow(omega/om_cp*(1.0 - M)*cos(theta), 2));
      
      //Steal Psi value from mu using Lyons \cite Lyons1974B Eq 12
      //Psi from "warm" theory, high dens, see mu assignment
      Psi = std::pow(om_pe_local/om_ce_local, 2) * (1.0+M)/M / my_mu.mu/my_mu.mu;

      g_x = my_spect->get_g_element(i);

      I_contrib = x * std::pow( (1.0 + x*x)*Psi, -1.5);
      I_contrib2 = 1.0 + 1.0/Psi * (om_sq_p_e - 0.5 * std::pow(omega/om_pe_local*(1.0-M), 2)/(( 1.0+ x*x)*(Psi - 1.0 + om_sq_p_e ) + 0.5*x*x ) );

      I_om += g_x * I_contrib * I_contrib2 * dx;
    }
  }
  return I_om;
}

//----------------------------------------------------------------

test_entity_levelone::test_entity_levelone(){
/** \brief Setup level-one testing
\todo Create some lighter weight test files!, perhaps SDF filter??*/
  name = "level-one derivation";
  
  test_contr = nullptr;
  my_reader = nullptr;
  
}
test_entity_levelone::~test_entity_levelone(){
/** \brief Teardown level-one testing*/
  delete test_contr;
  delete my_reader;
}

int test_entity_levelone::run(){
/** \brief Test entire level-1 data extraction
*
*Tests the full sequence of file reading, FFT and spectrum generation. Effectively a regression test as we check against previously generated files
*
Set runtime_flag "level_one" to perform full level-one testing
*
@return Error code
**/

  int err = TEST_PASSED;

  get_deck_constants(tests_src_dir + "test");
  share_consts();

  if(test_bed->runtime_flags.count("level_one") != 0){
    block_id = "ay";
    file_prefix = tests_src_dir + "l1/l1";
    space_in[0] = 0;
    space_in[1] = 1024;
    time_in[0] = 0;
    time_in[1] = 4;
    time_in[2] = 100;

    err|= setup();
    if(!test_bed->check_for_abort(err)) err|= basic_tests(1, -1, true);
    if(my_reader){
      delete my_reader;
      my_reader = nullptr;
    }
    if(test_contr){
      delete test_contr;
      test_contr = nullptr;
    }

    if(!test_bed->check_for_abort(err)){
      file_prefix = tests_src_dir + "2dtest/";
      block_id = "ey";
      space_in[0] = 0;
      space_in[1] = 1024;
      time_in[0] = 0;
      time_in[1] = 50;
      time_in[2] = 0;

      err|=setup();
      if(!test_bed->check_for_abort(err)) err |= basic_tests(2, 1, true);
      if(!test_bed->check_for_abort(err)) err |= basic_tests(2, -1, false, "_space", 2, 0.01f*my_const.omega_ce, 1.5f*my_const.omega_ce);
    }
  }else{
    test_bed->report_info("Skipping level-one tests, flag -level_one not set", 0);
  }
  return err;
}

int test_entity_levelone::setup(){
/** \brief Setup for "level one" extraction
*
*Creates required file reader and controller objects
@return Error code
*/

  int err = TEST_PASSED;
  my_reader = new reader(file_prefix, block_id);
  if(my_reader->current_block_is_accum()){
    n_tims = time_in[2];
  }else{
    n_tims = std::max((int) (time_in[1]-time_in[0]), 1);
  }

  size_t n_dims;
  std::vector<size_t> dims;
  int err2 = my_reader->read_dims(n_dims, dims);
  if(err2) err |= TEST_FATAL_ERR;
  
  test_contr = new controller(file_prefix);

  return err;
}

int test_entity_levelone::basic_tests(size_t n_dims_in, int flatten_on, bool has_freq, std::string outfile_tag, int total_fft, my_type band_min, my_type band_max){
/** \brief Basic tests of process to make level-1 data
*
* Reads proper data files, produces FFT, and derived spectrum, writes to file. Compares to reference version for regression. 
@param n_dims_in is spatial dimension expected for input file 
@param flatten_on Flattening dimension on raw data, negative for no flattening 
@param has_freq Output array contains frequency as last axis 
@param outfile_tag Tag for output file if wanted 
@param total_fft Total fft on this dimension, negative for no totalling 
@param band_min Minimum of band to total fft on
@param band_max Maximum of band to total fft on
@return Error code*/

  int err = TEST_PASSED;

  size_t n_dims;
  std::vector<size_t> dims;
  int err2 = my_reader->read_dims(n_dims, dims);
  if(err2){
    err |= TEST_FATAL_ERR;
    return err;
  }
  if(n_dims != n_dims_in){
    test_bed->report_info("Wrong file dimension", 1);
    return TEST_FATAL_ERR;
  }

  //Cut out section here
  int space_dim = space_in[1]-space_in[0];

  //Different size for space processing versus "normal"
  if(total_fft >= 0){
    dat = data_array(space_dim, dims[1], n_tims);
  }else{
    dat = data_array(space_dim, n_tims);
  }
  if(!dat.is_good()){
    my_error_print("Data array allocation failed.", mpi_info.rank);
    err |= TEST_ASSERT_FAIL;
    err |= TEST_FATAL_ERR;
  }
  if(flatten_on < 0){
    err2 = my_reader->read_data(dat, time_in, space_in);
  }else{
    err2 = my_reader->read_data(dat, time_in, space_in, flatten_on);
  }

  if(err2 == 1){
    return TEST_FATAL_ERR;
  }
  if(err2 == 2) n_tims = dat.get_dims(dat.get_dims() - 1);
  //Check if we had to truncate data array...

  dat_fft = data_array();
  dat_fft.clone_empty(dat);
  if(!dat_fft.is_good()){
    return TEST_FATAL_ERR;
  }
  dat.B_ref = get_ref_Bx(file_prefix, space_in, 1);
  //Use first file. 0th accumulator is perhaps broken
  err2 = fft_array(dat,dat_fft);

  if(total_fft >= 0){
    dat_fft = dat_fft.total(total_fft, band_min, band_max);
  }

  test_bed->report_info("FFT returned err_state " + mk_str(err2));

  test_contr->set_plasma_B0(dat.B_ref);
  test_contr->add_spectrum(space_dim, DEFAULT_N_ANG, true);
  test_contr->get_current_spectrum()->generate_spectrum(dat_fft, 10, FUNCTION_DELTA, DEFAULT_SPECTRUM_ANG_STDDEV);
  
  n_dims = dat_fft.get_dims();
  //Set cutout limits on FFT
  std::vector<my_type> lims;
  if(n_dims >=3){
    lims.push_back(-0.002);
    lims.push_back(0.002);
  }
  if(n_dims >=2){
    lims.push_back(-0.002);
    lims.push_back(0.002);
    if(has_freq){
      lims.push_back(-3.0*my_const.omega_ce);
      lims.push_back(3.0*my_const.omega_ce);
    }else{
      lims.push_back(-0.002);
      lims.push_back(0.002);
    }
  }
  
  //Dump files and then compare to reference files
  std::string filename, full_filename, time_str;
  time_str = mk_str(dat_fft.time[0], true)+"_"+mk_str(this->n_tims);
  std::string block = block_id;
  filename = "FFT_"+block +"_"+time_str+"_"+mk_str(dat_fft.space[0])+"_"+mk_str(dat_fft.space[1]) + outfile_tag+".dat";
  full_filename = tests_tmp_dir + filename;
  std::fstream file;
  file.open(full_filename.c_str(),std::ios::out|std::ios::binary);
  if(file.is_open()){
    dat_fft.write_section_to_file(file, lims);
    if(err2){
      test_bed->report_info("File writing failed");
      err |=TEST_ASSERT_FAIL;
    }
    
  }else{
    err |= TEST_ASSERT_FAIL;
  
  }
  file.close();
  test_bed->report_info("FFT section output in "+full_filename, 1);
  data_array previous_fft(file_prefix + "tests/"+ filename+".ref");
  dat_fft = data_array(full_filename);
  if(previous_fft != dat_fft){
    test_bed->report_info("New FFT does not match reference");
    err |= TEST_WRONG_RESULT;
  }

  full_filename = append_into_string(full_filename, "_spectrum");
  filename = append_into_string(filename, "_spectrum");
  file.open(full_filename.c_str(),std::ios::out|std::ios::binary);
test_contr->get_current_spectrum()->wave_id = 0;
  if(file.is_open()){
    test_contr->get_current_spectrum()->write_to_file(file);
    if(err2){
      test_bed->report_info("File writing failed");
      err |=TEST_ASSERT_FAIL;
    }
    
  }else{
    err |=TEST_ASSERT_FAIL;
  
  }
  file.close();
  test_bed->report_info("Spectrum output in "+full_filename, 1);

  test_contr->add_spectrum(file_prefix + "tests/" + filename + ".ref");
  if(test_contr->get_current_spectrum() != nullptr && test_contr->get_spectrum_by_num(1) != nullptr && *(test_contr->get_current_spectrum()) != *(test_contr->get_spectrum_by_num(1))){
    test_bed->report_info("New spectrum does not match reference");
    err |= TEST_WRONG_RESULT;
  }

  return err;

}

//----------------------------------------------------------------

test_entity_d::test_entity_d(){
/**  \brief Setup tests for dffusion calculation
*/
  name = "D checks";
  file_prefix = tests_src_dir + "d_test";
  test_contr = nullptr;
  partial_disable_for_valgrind = false;

}
test_entity_d::~test_entity_d(){
/** \brief Teardown tests for diffusion calcs*/
  if(test_contr) delete test_contr;
}

int test_entity_d::run(){
/** \brief Test D calculation
*
*A full useful calculation of a test coefficient takes quite a while, so here we just test that the calculation proceeds and there are no obvious problems.
*Set runtime_flag "full_d" to perform a full sample D calculation
@return Error code
*/

  int err = TEST_PASSED;
#ifdef RUNNING_ON_VALGRIND
  if(RUNNING_ON_VALGRIND){
    my_error_print("Detected valgrind.\n There is an extended precision bug so boost::math bessel functions don't work correctly.\n        Some D tests will not be performed", mpi_info.rank);
    partial_disable_for_valgrind = true;
  }
//Valgrind doesn't support extended precision, so the boost::math definition of a "tiny" to protect against divide by zero is equal 0 resulting in NaN. This is a known "bug" so we just skip certain D tests if valgrind is in use, as they will fail
#endif
  
  deck_constants const_tmp = my_const;
  if(mpi_info.rank == 0) get_deck_constants(file_prefix);
  share_consts();

  test_contr = new controller(file_prefix);
  err |= basic_tests();
  if(test_bed->runtime_flags.count("full_d") > 0){
    err |= full_D_tests();
  }else{
    test_bed->report_info("Skipping full_D tests. Invoke with flag -full_d", 0);
  }

  my_const = const_tmp;
  share_consts();

  return err;
}

int test_entity_d::basic_tests(){
/** \brief Simple tests of D
*
* Does some simple checks that D calc proceeds and some basic things are true. Also does an IO test
@return Error code
*/
  int err = TEST_PASSED;
  
  test_bed->report_info("Reading spectrum", mpi_info.rank);
  //Now dump to file and read back in and compare
  bool err2 = test_contr->add_spectrum(tests_tmp_dir + "spect_out.dat");
  if(!err2){
    test_bed->report_info("Calculating test D", mpi_info.rank);
    test_contr->add_d(5, 5);
    if(!test_contr->get_current_d()){
      //Can't go on, have no D!!!
      err |= TEST_FATAL_ERR;
      return err;
    }
    d_report report = test_contr->get_current_d()->calculate();
    if(report.error){
      test_bed->report_info("Error calculating D", mpi_info.rank);
      err |= TEST_ASSERT_FAIL;
    }

  }else{
    err |= TEST_ASSERT_FAIL;
  }

  test_bed->report_info("Testing D IO");
  //Modify some fields to non-default
  test_contr->get_current_d()->wave_id = 3;
  test_contr->get_current_d()->tag = "IO test, long tag";
  
  //Write the current D to file
  std::fstream outfile;
  std::string filename = tests_tmp_dir + "d_dump.txt";
  outfile.open(filename.c_str(), std::ios::binary|std::ios::out|std::ios::in|std::ios::trunc);
  test_contr->get_current_d()->write_to_file(outfile);
  outfile.close();
  //Create new empty D
  err2 = test_contr->add_spectrum(tests_tmp_dir + "spect_out.dat");
  test_contr->add_d(5, 5);
  //Read from file
  outfile.open(filename.c_str(), std::ios::binary|std::ios::in);
  test_contr->get_current_d()->read_from_file(outfile);
  //Check equality

  bool D_is_eq = true;//Flag for equality of written and read D
  bool D_is_bad = false;//Flag for if written D contains NaN elements
  bool D_bad_element = false;
  for(int i=0; i< 5; i++){
    for(int j=0; j< 5; j++){
      if(test_contr->get_d_by_num(1)->get_element(i, j) != test_contr->get_d_by_num(1)->get_element(i, j)){
         D_bad_element = true;//Element is NaN in written array
         D_is_bad = true;
      }
      if(!D_bad_element && test_contr->get_d_by_num(1)->get_element(i, j) != test_contr->get_current_d()->get_element(i, j)) D_is_eq = false;
        //We wrote a NaN element so will never be equal
      D_bad_element = false;
    }
  }
  if(test_contr->get_d_by_num(1)->wave_id != test_contr->get_current_d()->wave_id) D_is_eq = false;
  //Note tag gets truncated to ten chars in writing
  if(test_contr->get_d_by_num(1)->tag.substr(0, 9) != test_contr->get_current_d()->tag.substr(0, 9)) D_is_eq = false;

  if(D_is_bad && !partial_disable_for_valgrind){
    //Due to valgrind issue described in ::run(), D might be NaN but this is not our fault...
    test_bed->report_info("D contains NaN!");
    err |= TEST_WRONG_RESULT;
  }
  
  if(!D_is_eq){
    test_bed->report_info("Error reading D from file");
    err |= TEST_USERDEF_ERR1;
  }
  return err;
}

int test_entity_d::full_D_tests(){
/** \brief Calculate a test D
*
* Calculate a real diffusion coefficient for useful parameters, and check features against sample data
*Set runtime_flag "n" to a value to calculate a single n
*Set runtime_flag "d_pp" to calculate D_pp rather than alpha_alpha by default

@return Error code
*/

  int err = TEST_PASSED;
  bool single_n = false;
  int n = 0;

  bool err2 = test_contr->add_spectrum(file_prefix + "spectrum.dat");
  if(!err2){
    test_bed->report_info("Calculating full D... This may take a (very) long time!\nEnsure optimisation is on during compile!", mpi_info.rank);
    test_contr->add_d(100, 100);
    
    if(test_bed->runtime_flags.count("n") != 0){
      single_n = true;
      n  = test_bed->runtime_flags["n"];
      test_contr->get_current_d()->set_single_n(n);
    }else{
      test_contr->get_current_d()->set_max_n(5);
    }
    D_type_spec D_type = D_type_spec::alpha_alpha;
    if(test_bed->runtime_flags.count("d_pp") != 0){
      D_type = D_type_spec::p_p;
    }
    d_report report = test_contr->get_current_d()->calculate(D_type);
    if(report.error){
      test_bed->report_info("Error calculating full D", mpi_info.rank);
      err |= TEST_ASSERT_FAIL;
    }
    test_bed->report_info("Average resonant number "+mk_str(report.n_av), 2);//Num resonances is 2* this +1 
    if(!report.single_n) test_bed->report_info("Max resonant number "+mk_str(report.n_max), 2);
    
    test_bed->report_info("Writing test file", mpi_info.rank);
    std::fstream file;
    file.open("test_d.dat", std::ios::binary|std::ios::trunc|std::ios::out|std::ios::in);
    if(file) test_contr->get_current_d()->write_to_file(file);
    file.close();

    data_array padie_data = read_padie_data(single_n, n, (D_type == D_type_spec::p_p));
    //Now we should read a file containing sample D info and compare some features. Perhaps the Landau peak?
    //Check axes match. Remember D has first velocity, padie is just 1-d of angle
    if(padie_data.is_good()){
      //Find energy bin to use. Padie is at 1 MeV, so 0.941c
      size_t v_ind = test_contr->get_current_d()->get_axis_index_from_value(0, 0.941*v0);
      calc_type local_diff, running_diff = 0.0, running_total = 0.0;
      //Calculate the ~L^2 difference
      for(size_t i = 0; i < padie_data.get_dims(0); i++){
        local_diff = std::pow(padie_data.get_element(i) - test_contr->get_current_d()->get_element(v_ind, i), 2);
        running_diff += local_diff;
        running_total += std::pow(padie_data.get_element(i), 2);
      }
      test_bed->report_info("Calculated D has L2 mismatch of " + mk_str((int)(std::sqrt(running_diff/running_total)*100 ))+'%');
      //Very relaxed error condition
      if(std::sqrt(running_diff/running_total)*100 > 50) err |= TEST_WRONG_RESULT;
    }else{
      err |= TEST_ASSERT_FAIL;
    }
    
  }else{
    err |= TEST_ASSERT_FAIL;
  }


  return err;
}

data_array test_entity_d::read_padie_data(bool single_n, int n, bool use_d_pp){
/** \brief Read reference data
*
* Read data from reference PADIE output files into a new data array
@param single_n Read data for a single resonance
@param n The resonant n to read if single_n set
@param use_d_pp Return d_pp data rather than d_alpha alpha
@return New data array containing reference data
*/
  std::string padie_dir = tests_src_dir + "PadieTestData/";
  std::string filestart = "Data";
  std::string fileend = "", filename;

  if(single_n){
    fileend = "_";
    fileend += mk_str(n);
  }
  filename = padie_dir + filestart + fileend;
  if(use_d_pp) filename += "_p.dat";
  else filename += ".dat";
  test_bed->report_info("Reading " + filename);
  data_array data = data_array(filename);
  
  return data;
}

//----------------------------------------------------------------

test_entity_bounce::test_entity_bounce(){
/** \brief Setup tests for bounce averaging*/
  name = "bounce averaging";
  file_prefix = tests_src_dir;
  //Have to set-up deck constants before creating plasma
  get_deck_constants(tests_src_dir + "test");
  share_consts();

  test_contr = new controller(file_prefix);
}
test_entity_bounce::~test_entity_bounce(){
/** \brief Teardown tests for bounce averaging*/
  if(test_contr) delete test_contr;
}

int test_entity_bounce::run(){
/** \brief Test bounce averaging
*
* Tests the 4 forms of bounce averaging, simple line-length totalling and the p_p, alpha_p and alpha_alpha forms
@return Error code
*/

  int err = TEST_PASSED;

  //First we do a simple test by rigging up some space-blocked 1x1 D's containing value 1, This should return the length of the line which we calc using Eq 1.01 in Schulz/Lanzerotti and compare

  bounce_av_data bounce_dat;
  bounce_dat.max_latitude = 90.0;
  bounce_dat.L_shell = 4.5;

  my_type line_length;
  line_length = 2.7603/2.0 * bounce_dat.L_shell * R_E;
  test_bed->report_info("Expecting line length of "+mk_str(line_length/1000.0)+" km");

  for(size_t n=0; n<32; n++){
    test_contr->add_spectrum(1, 1, true);
    test_contr->add_d(1, 1);
    test_contr->get_current_d()->set_element((size_t)0,(size_t)0, 1.0);
  }
  test_contr->bounce_average(bounce_dat);
  //These should match to within say 0.1% for a reasonable result
  test_bed->report_info("Got line length "+mk_str(test_contr->get_special_d()->get_element((size_t)0,(size_t)0)/1000.0)+" km");

  if(std::abs(test_contr->get_special_d()->get_element((size_t)0,(size_t)0)/line_length - 1.0) > 0.001){
    err |= TEST_WRONG_RESULT;
    test_bed->report_info("Line length mismatch of "+mk_str(std::abs(test_contr->get_special_d()->get_element((size_t)0,(size_t)0)/line_length - 1.0)*100.0, false) +" %");
  }
  
  //Now we check over some of the helpers, mirror latitude, bounce period etc using a sample pitch angle and latitude
  my_type test_alpha_eq = 10.0*pi/180.0, expected_mirror_lat = 52, expected_bounce_period = 1.2, expected_alpha_at_lat_20 = 13, tmp;
  tmp = solve_mirror_latitude(test_alpha_eq)*180/pi;
  if(std::abs(tmp - expected_mirror_lat) > 1.0){
    //Mirror lat to nearest degree
    err |= TEST_WRONG_RESULT;
    test_bed->report_info("Wrong mirror latitude, value "+mk_str(tmp)+" expected "+mk_str(expected_mirror_lat), 1);
  }
  tmp = bounce_period_approx(test_alpha_eq);
  if(std::abs(tmp - expected_bounce_period) > 0.01){
    //Bounce period factor to 3 dp
    err |= TEST_WRONG_RESULT;
    test_bed->report_info("Wrong bounce period, value "+mk_str(tmp)+" expected "+mk_str(expected_bounce_period), 1);
  }
  tmp = alpha_from_alpha_eq(test_alpha_eq, 20.0*pi/180.0)*180/pi;
  if(std::abs(tmp - expected_alpha_at_lat_20) > 1.0){
    //Pitch angle at lat 20 deg. to nearest degree
    err |= TEST_WRONG_RESULT;
    test_bed->report_info("Wrong pitch angle, value "+mk_str(tmp)+" at latitude 20 degrees, expected "+mk_str(expected_alpha_at_lat_20), 1);
  }

  if(err == TEST_PASSED) test_bed->report_info("Bounce helpers OK", 2);
  
  bounce_dat.type = p_p;
  err |= bounce_cases(bounce_dat);
  bounce_dat.type = alpha_alpha;
  err |= bounce_cases(bounce_dat);
  bounce_dat.type = alpha_p;
  err |= bounce_cases(bounce_dat);

  return err;

}

int test_entity_bounce::bounce_cases(bounce_av_data bounce_dat){
/** \brief Test bounce averaging
*
* Tests the p_p, alpha_p and alpha_alpha bounce averages by setting up a D so that the bounce average contains no alpha dependency, only lambda and alpha_eq and eval. This gives integrals that can be analytically integrated, which we do and check, allowing a "moderate" discrepancy given the small number of space blocks
@param bounce_dat Bounce averaging data such as type and L-shell
@return Error code
*/
  
  int err = TEST_PASSED;
  test_contr->clear_all();
  size_t d_sz =  150;
  my_type val;
  for(size_t n=0; n<32; n++){
    test_contr->add_spectrum(1, 1, true);
    test_contr->add_d(d_sz, d_sz);
    //Do axes
    for(size_t i=0; i< d_sz; i++){
      for(size_t j=0; j< d_sz; j++){
        switch(bounce_dat.type){
          case plain:
            val = 0.0;
            break;
          case p_p:
            val = cos(test_contr->get_current_d()->get_axis_element_ang(j));
            break;
          case p_alpha:
          case alpha_p:
            val = 1.0/sin(test_contr->get_current_d()->get_axis_element_ang(j));
            //At zero we'll use any finite value because the numerator will be zero
            if(test_contr->get_current_d()->get_axis_element_ang(j) == 0.0) val = 1.0;
            break;
          case alpha_alpha:
            val = 1.0/cos(test_contr->get_current_d()->get_axis_element_ang(j));
            break;
        }
        test_contr->get_current_d()->set_element(i,j, val);
      }
    }
  }
  test_contr->bounce_average(bounce_dat);
  //Check the value at random p and all angles
  my_type lambda_m, current_val, expected_val, alpha_eq;
  std::string err_string;
  size_t i = 10;
  //Cheating and bumping down because at large initial pitch angles we get it wrong a bit
  /** \todo Can we fix this at large angle, or small angle for sin?*/
  for(size_t j=1; j< d_sz; j++){
    alpha_eq = test_contr->get_current_d()->get_axis_element_ang(j);
    lambda_m = solve_mirror_latitude(alpha_eq);
    current_val = test_contr->get_special_d()->get_element(i,j);
    
    switch(bounce_dat.type){
      case plain:
        return TEST_ASSERT_FAIL;//Not meant to try this here
      case p_p:
        expected_val =(3.0*sin(lambda_m)*std::sqrt(3.0*std::pow(sin(lambda_m), 2) + 1.0) + std::sqrt(3.0)*asinh(std::sqrt(3.0)*sin(lambda_m)))/6.0/bounce_period_approx(alpha_eq);
        err_string = "p_p";
        if(j > d_sz - 12) err |= TEST_REMOVE_ERR;
        break;
      case p_alpha:
      case alpha_p:
        expected_val = (1225.0*sin(lambda_m) + 245.0*sin(3.0*lambda_m) + 49.0*sin(5.0*lambda_m) + 5.0*sin(7.0*lambda_m))/2240.0/cos(alpha_eq)/sin(alpha_eq)/bounce_period_approx(alpha_eq);
        err_string = "alpha_p";
        if(j < 12) err |= TEST_REMOVE_ERR;
        break;
      case alpha_alpha:
        expected_val = (1225.0*sin(lambda_m) + 245.0*sin(3.0*lambda_m) + 49.0*sin(5.0*lambda_m) + 5.0*sin(7.0*lambda_m))/2240.0/std::pow(cos(alpha_eq), 2)/bounce_period_approx(alpha_eq);
        err_string = "alpha_alpha";
        if(j > d_sz - 12) err |= TEST_REMOVE_ERR;
        break;
    }
    if(std::abs(current_val/expected_val -1.0) > 0.05){
      //Moderate discrep allowed here
      err |= TEST_WRONG_RESULT;
      if((err & TEST_REMOVE_ERR) != TEST_REMOVE_ERR) test_bed->report_info("Erroneous bounce average in "+err_string+" mismatch "+mk_str((int)(std::abs(current_val/expected_val -1.0)*100))+'%', 2);
    }
  }
  return err;
}

test_entity_nonthermal::test_entity_nonthermal(){
/** \brief Setup non-thermal tests*/
  name = "nonthermal";
}

int test_entity_nonthermal::run(){
/** \brief Check nonthermal module
*
* Check function binding, df/dp values, lookup table
@return Error code
*/
  int err = TEST_PASSED;
  
  //Setup single Max non-thermal
  non_thermal * my_elec = new non_thermal(tests_src_dir + "test");

  calc_type a_par = std::sqrt(2.0)*my_elec->get_v_par() / std::sqrt(1.0 - std::pow(my_elec->get_v_par()/v0, 2));
  calc_type a_perp_sq = 2.0*std::pow(my_elec->get_v_perp(), 2)/(1.0 - 2.0*std::pow(my_elec->get_v_perp()/v0, 2));

  calc_type norm = (my_elec->get_total_dens() - 1.0)/pi/std::sqrt(pi)/a_par/a_perp_sq;

  calc_type df_tmp, f_tmp = my_elec->f_p(0, 0);
  
  if(std::abs(f_tmp / norm -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  const size_t n_trials = 5;
  calc_type inc = 0.1, p_par, p_perp;

  //Test f against expected values
  for(size_t i=0; i<n_trials*n_trials; i++){
    p_perp = inc* (i % n_trials);
    p_par = inc* (int)(i/n_trials);
    f_tmp = my_elec->f_p(p_par*a_par, p_perp*std::sqrt(a_perp_sq));
    if(std::abs(f_tmp /( norm*std::exp(-p_par*p_par)*std::exp(-p_perp*p_perp)) -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  }
  //Test a couple of extreme cases
  p_perp = 1.5;
  f_tmp = my_elec->f_p(p_par*a_par, p_perp*std::sqrt(a_perp_sq));
  if(std::abs(f_tmp /( norm*std::exp(-p_par*p_par)*std::exp(-p_perp*p_perp)) -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  p_par = 2.5;
  f_tmp = my_elec->f_p(p_par*a_par, p_perp*std::sqrt(a_perp_sq));
  if(std::abs(f_tmp /( norm*std::exp(-p_par*p_par)*std::exp(-p_perp*p_perp)) -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

  //Test f against expected values
  for(size_t i=0; i<n_trials; i++){
    p_perp = inc*i + inc;
    p_par = inc*i + inc;
    //close to 0 deriv -> 0

    f_tmp = my_elec->d_f_p(p_par*a_par, 0, 1);
    if(std::abs(f_tmp /((-2.0)*p_par/a_par*my_elec->f_p(p_par*a_par, 0)) -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;

    f_tmp = my_elec->d_f_p(0, p_perp*std::sqrt(a_perp_sq), 0);
    if(std::abs(f_tmp /((-2.0)*p_perp/std::sqrt(a_perp_sq)*my_elec->f_p(0, p_perp*std::sqrt(a_perp_sq))) -1.0) > GEN_PRECISION) err |=TEST_WRONG_RESULT;
  }

  //This next is deep wings so we relax constraint so we're not trying to resolve 1e-60
  df_tmp = my_elec->d_f_p(0, 10.0*std::sqrt(a_perp_sq), 0);
  if(std::abs(df_tmp /((-2.0)*10.0/std::sqrt(a_perp_sq)*my_elec->f_p(0, 10.0*std::sqrt(a_perp_sq))) -1.0) > 2e-6) err |=TEST_WRONG_RESULT;
  
  if(err == TEST_PASSED) test_bed->report_info("Nonthermal values OK",1);
  
  err |= test_lookup();
  return err;
}

int test_entity_nonthermal::test_lookup(){
/** \brief Check lookup table
*
*Check known lookup function against known analytic results. Data generated by generate_lookup in generate_test_data.pro
@return Error code
*/
  int err = TEST_PASSED;
  
  non_thermal * my_elec_an = new non_thermal(tests_src_dir + "an");
  non_thermal * my_elec_lookup = new non_thermal(tests_src_dir + "lookup");

  calc_type f_an, f_lookup, p_par, p_perp;
  const size_t n_tests = 100;

  for(size_t i=0; i< n_tests; i++){
    for(size_t j=0; j< n_tests; j++){
      p_par = (float) i *v0/2.0/(float) n_tests;
      p_perp = (float) j *v0/2.0/(float) n_tests;
      f_an = my_elec_an->f_p(p_par, p_perp);
      f_lookup = my_elec_lookup->f_p(p_par, p_perp);
      if(std::abs(f_lookup/f_an - 1.0) > LOW_PRECISION){
        err |= TEST_WRONG_RESULT;
      }
    }
  }
  if(err == TEST_PASSED) test_bed->report_info("Lookup OK", 1);
  return err;
}

#endif