alpha_calc.c

Go to the documentation of this file.

/** @file */
/*
** alpha_calc.c
** 
** made by (Johnny Q. Hacker)
** Login   <solarion@borkborkbork>

Copyright (C) 2009 Joseph Pingenot

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Affero General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Affero General Public License for more details.

You should have received a copy of the GNU Affero General Public License
along with this program.  If not, see
<http://www.gnu.org/licenses/>.

Equation is based off of JApplPhys_83_4324

** Started on  Fri Mar  6 11:22:01 2009 Johnny Q. Hacker
** Last update Wed Oct  7 15:15:00 2009 Johnny Q. Hacker
*/

#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <math.h>
#include <libcommon/gimme-alpha_structs.h>
#include <libpostproccalx/alpha_calc.h>
#include <libmodeling/wavefunction.h>
#include <libcommon/gcc_hints.h>

        /*Argument to the GSL integration workspace allocation calls*/
#define MAX_GSL_INTEGRATION_INTERVALS 4096


extern const double hbar;
extern const double q_e;
extern const double oo_q_e;
extern const double m_e_0;
extern const double oo_m_e_0;

/*eV, for InSb*/
const double Ep_InSb=23.3;
/*eV, for AlSb*/
const double Ep_AlSb=18.7;
/*eV, for InSb*/
const double Gamma8_edge_InSb=23.3;
/*eV, for AlSb*/
const double Gamma8_edge_AlSb=18.7;
/*eV, for InSb*/
const double Gamma7_edge_InSb=23.3;
/*eV, for AlSb*/
const double Gamma7_edge_AlSb=18.7;

/*For the InGaAs; eV*/
const double Ep_InAs=21.5;
const double Ep_GaAs=28.8;

/*x is % InAs*/
__const double InGaAs_Ep(double x) {
  return Ep_GaAs + (Ep_InAs - Ep_GaAs)*x;
}

___const double AlInSb_gamma8_bandedge(float x) {
  return Gamma8_edge_AlSb + (Gamma8_edge_InSb - Gamma8_edge_AlSb)*x;
}

___const double AlInSb_gamma7_bandedge(float x) {
  return Gamma7_edge_AlSb + (Gamma7_edge_InSb - Gamma7_edge_AlSb)*x;
}

___const double AlInSb_Ep(float x) {
  return Ep_AlSb + (Ep_InSb - Ep_AlSb)*x;
}

___const double complex_magnitude(double re, double im) {
  return re*re + im*im;
}

/*Information for the gsl integration callback.*/
struct expectation_value_params {
  gsl_interp_accel *re_accel;
  gsl_interp_accel *im_accel;
  gsl_interp_accel *func_accel;
  gsl_spline *re;
  gsl_spline *im;
  gsl_spline *func;
};
/*This part is simple.  We just want <psi|func|psi>, where func is
  ONLY a function of x
  *NOTE: this is provided as a callback to gsl.
  *NOTE: there are 3 flavours; one for each derivative.
  */
double expectation_value_integrand0(double x, void* params) {
  return complex_magnitude(
         gsl_spline_eval(((struct expectation_value_params*)params)->re, x, ((struct expectation_value_params*)params)->re_accel),
         gsl_spline_eval(((struct expectation_value_params*)params)->im, x, ((struct expectation_value_params*)params)->im_accel)
         )*gsl_spline_eval(((struct expectation_value_params*)params)->func, x, ((struct expectation_value_params*)params)->func_accel);
}
double expectation_value_integrand1(double x, void* params) {
  return complex_magnitude(
         gsl_spline_eval(((struct expectation_value_params*)params)->re, x, ((struct expectation_value_params*)params)->re_accel),
         gsl_spline_eval(((struct expectation_value_params*)params)->im, x, ((struct expectation_value_params*)params)->im_accel)
         )*gsl_spline_eval_deriv(((struct expectation_value_params*)params)->func, x, ((struct expectation_value_params*)params)->func_accel);
}
double expectation_value_integrand2(double x, void* params) {
  return complex_magnitude(
         gsl_spline_eval(((struct expectation_value_params*)params)->re, x, ((struct expectation_value_params*)params)->re_accel),
         gsl_spline_eval(((struct expectation_value_params*)params)->im, x, ((struct expectation_value_params*)params)->im_accel)
         )*gsl_spline_eval_deriv2(((struct expectation_value_params*)params)->func, x, ((struct expectation_value_params*)params)->func_accel);
}

/*!\brief calculates the expectation value of func in state wf.
 *ARGS:
 * \param function: array of N doubles containing the function.
 * \paramx: array of N doubles containing the x points at which function was sampled
 * \param N: the numnber of points in discrete function function.
 * \param wf: pointer to an mstar_wavefunction_1d structure containing
 * \param derivative: which derivative, if any, to use in evaluating the
 *   function, e.g. derivative=1 -> find the expectation value of the
 *   first derivative of function in state wf.
 * \param abserr largest acceptable absolute integration error (GSL)
 * \param relerr largest acceptable relative integration error (GSL)
 * \param err pointer to a double to store the value of the (GSL)
 *   estimated integration error.
 *RETURNS:
 * The expectation value of the derivative-th derivative of function
 *in state wf.
 *NOTE: function must ONLY be a function of x!
 *NOTE: only up to second derivative is supported.
 */
double average_discrete_function(const double *function, const double* x, int N, const struct mstar_wavefunction_1d *wf, short int derivative, double abserr, double relerr, double* err) {
  gsl_interp_accel *func_accel = gsl_interp_accel_alloc();
  gsl_spline *func = gsl_spline_alloc(gsl_interp_akima, N);
  gsl_spline_init(func, x, function, N);
  /*Wavefunction*/
  double *repacked_x;
  double *repacked_re;
  double *repacked_im;
  mstar_wavefunction_1d_repack(wf, &repacked_x, &repacked_re, &repacked_im, 1);
  gsl_interp_accel *re_accel = gsl_interp_accel_alloc();
  gsl_spline *re = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(re, repacked_x, repacked_re, wf->N);
  gsl_interp_accel *im_accel = gsl_interp_accel_alloc();
  gsl_spline *im = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(im, repacked_x, repacked_im, wf->N);
  struct expectation_value_params expectation_value_params = {
    .re_accel = re_accel,
    .im_accel = im_accel,
    .func_accel = func_accel,
    .re = re,
    .im = im,
    .func = func
  };
  gsl_function gsl_func = {
    .function = &expectation_value_integrand0,
    .params = &expectation_value_params
  };
  double e_val;
  double status=0;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(MAX_GSL_INTEGRATION_INTERVALS);
  switch(derivative) {
  case 0:
    gsl_func.function = &expectation_value_integrand0;
    break;
  case 1:
    gsl_func.function = &expectation_value_integrand1;
    break;
  case 2:
    gsl_func.function = &expectation_value_integrand2;
    break;
  default:
    return nan("char-sequence");
  }
  /*Here's where we do the integration.*/
  status = gsl_integration_qag(&gsl_func, repacked_x[0], repacked_x[wf->N-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w, &e_val, err);
  gsl_spline_free(re);
  gsl_spline_free(im);
  gsl_spline_free(func);
  gsl_interp_accel_free(re_accel);
  gsl_interp_accel_free(im_accel);
  gsl_interp_accel_free(func_accel);
  free(repacked_x);
  free(repacked_re);
  free(repacked_im);
  return e_val;
}

/*! Averages n_arrays of array_size doubles.
 *
 *\param arrays  pointer to an array of arrays of doubles
 *\param n_arrays  number of arrays given
 *\param array_size the size of the arrays stored in arrays.
 */
double* average_arrays(double *arrays[], int n_arrays, int array_size) {
  if(unlikely(arrays == NULL)) return NULL;
  double *avg = malloc(sizeof(double)*array_size);
  if(unlikely(avg == NULL)) return NULL;
  if(unlikely(n_arrays < 1)) {
    free(avg);
    return NULL;
  }
  register int i;
  for(i=0; i<array_size; i++) {
    avg[i] = (arrays[0])[i];
    //fprintf(stderr, "(arrays[0])[%d]=%g\tavg[%d]=%g\n", i, (arrays[0])[i], i, avg[i]);
  }
  for(i=1; i<n_arrays; i++) {
    for(int j=0; j<array_size; j++) {
      avg[j] += (arrays[i])[j];
      //fprintf(stderr, "(arrays[%d])[%d]=%g\tavg[%d]=%g\n", i, j, (arrays[i])[j], i, avg[j]);
    }
  }
  for(i=0; i<array_size; i++) {
    avg[i] /= n_arrays;
  }
  return avg;
}

/*!Finds the points to the left and right of the specified interface
  \param iface  interface point
  \param prev   pointer to an int to store the index of the previous point
  \param next   Same, but the next point
  \param x      array of x values, assumed to be sorted from smallest to largest
  \param N      number of points in the x array
  \return pointer to a string that contains an error message (NULL if succeeded)
*/
char* find_left_right_of_iface(double iface, int* prev, int* next, double *x, int N) {
  register int i;
  for(i=0; i<N; i++) {
    if(x[i] >= iface) {
      *prev = i-1;
      if(x[i] == iface) {
  //fprintf(stderr, "=(p=%i, n=%i, %g==%g [%g < %g])\n", *prev, *next, x[i], iface, x[*prev], iface);
  *next = i+1;
      }else{
  //fprintf(stderr, ">(p=%i, n=%i, %g>%g [%g < %g])\n", *prev, *next, x[i], iface, x[*prev], iface);
  *next = i;
      }
      break;
    }
    //fprintf(stderr, ".");
  }
  if(unlikely((i == N) || (*prev < 0) || (*next >= N))) {
    //fprintf(stderr, "E\n");
    if(i == N) return "specified interface not in structure";
    if(*prev < 0) return "specified interface is the start of the structure";
    if(*next >= N) return "specified interface is the end of the structure";
  }
  return NULL;
}


double alpha_interface_term_mstar_nostrain(const struct mstar_wavefunction_1d *wf, const struct band *hh_band, const struct band *lh_band, const struct band *so_band, const struct band *cb_band, const struct potential *potential, GList* ifs, double Ep) {
  double alpha = 0;
  int i=0;
  /*HH*/
  gsl_interp_accel *hh_acc = gsl_interp_accel_alloc();
  gsl_spline *hh = gsl_spline_alloc(gsl_interp_akima, hh_band->N);
  gsl_spline_init(hh, hh_band->x, hh_band->bandedge, hh_band->N);
  /*LH*/
  gsl_interp_accel *lh_acc = gsl_interp_accel_alloc();
  gsl_spline *lh = gsl_spline_alloc(gsl_interp_akima, lh_band->N);
  gsl_spline_init(lh, lh_band->x, lh_band->bandedge, lh_band->N);
  /*SO*/
  gsl_interp_accel *so_acc = gsl_interp_accel_alloc();
  gsl_spline *so = gsl_spline_alloc(gsl_interp_akima, so_band->N);
  gsl_spline_init(so, so_band->x, so_band->bandedge, so_band->N);
  /*CB*/
  gsl_interp_accel *cb_acc = gsl_interp_accel_alloc();
  gsl_spline *cb = gsl_spline_alloc(gsl_interp_akima, cb_band->N);
  gsl_spline_init(cb, cb_band->x, cb_band->bandedge, cb_band->N);
  /*Potential*/
  gsl_interp_accel *pot_acc = gsl_interp_accel_alloc();
  gsl_spline *pot = gsl_spline_alloc(gsl_interp_akima, potential->N);
  gsl_spline_init(pot, potential->x, potential->potential, potential->N);
  /*Wavefunction*/
  double *repacked_x;
  double *repacked_re;
  double *repacked_im;
  //for(unsigned i=0; i<wf->N; i++) fprintf(stderr, "wf.x=%g wf.re=%g\n", wf->points[i].x, wf->points[i].re);
  mstar_wavefunction_1d_repack(wf, &repacked_x, &repacked_re, &repacked_im, 1);
  //for(unsigned i=0; i<wf->N; i++) fprintf(stderr, "->repacked_x=%g repacked_re=%g\n", repacked_x[i], repacked_re[i]);
  gsl_interp_accel *re_acc = gsl_interp_accel_alloc();
  gsl_spline *re = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(re, repacked_x, repacked_re, wf->N);
  gsl_interp_accel *im_acc = gsl_interp_accel_alloc();
  gsl_spline *im = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(im, repacked_x, repacked_im, wf->N);
  /*Evaluate at the interfaces*/
  double *arrays[] = {hh_band->bandedge, lh_band->bandedge};
  double *avg = average_arrays(arrays, 2, hh_band->N);
  int prev=0, next=0;
  char* err;
  /*We don't need to keep track of the sign of the interface; that's
   *preserved from the definition of A and B; we just blindly use prev
   *and next; see epsilon_7^A and friends below for why*/
  double avg_p;
  double avg_n;
  double so_p;
  double so_n;
  double pot_p;
  double pot_n;
  /*the 0-point of energies reported in the printf below*/
  int zero_index_prev=0, zero_index_next=0;
  double relative_to;
  if(unlikely((err = find_left_right_of_iface(0, &zero_index_prev, &zero_index_next, hh_band->x, hh_band->N)) != NULL)) {
    fprintf(stderr, "WARNING: diagnostic print band structure is going to be relative to 0 (E_Fermi) (err=%s)\n", err);
    relative_to=0;
  }else{
    /*Average prev and next if there's no actual point at 0.*/
    if((zero_index_prev + 1) == zero_index_next) {
      relative_to = (((cb_band->bandedge)[zero_index_prev+1]) + ((cb_band->bandedge)[zero_index_prev+1]))*0.5;
    }else{
      relative_to = (cb_band->bandedge)[zero_index_prev+1];
    }
  }
  for(GList *i=g_list_first(ifs); i!=NULL; i=g_list_next(i)) {
    if(unlikely((err = find_left_right_of_iface(*((double*)(i->data)), &prev, &next, hh_band->x, hh_band->N)) != NULL)) {
      fprintf(stderr, "ERROR when getting prev/next indices for interface at %g: %s\n", *((double*)(i->data)), err);
      return -0;
    }
    /*
      epsilon_7^A is so_band->bandedge)[next]
      epsilon_7^B is (so_band->bandedge)[prev]
      epsilon_8^A is avg[prev]
      epsilon_8^B is avg[next]
      delta_7 = epsilon_7^B - epsilon_7^A 
      = (so_band->bandedge)[next] - (so_band->bandedge)[prev]
      delta_8 = epsilon_8^B - epsilon_8^A 
      = avg[next] - avg[prev]
      NOTE that the above have the potential still in them!
    */
    avg_p = avg[prev];
    avg_n = avg[next];
    so_p = (so_band->bandedge)[prev];
    so_n = (so_band->bandedge)[next];
    pot_p = (potential->potential)[prev];
    pot_n = (potential->potential)[next];
    //fprintf(stderr, "%g: d_avg=%g d_so=%g avg_n=%g avg_p=%g so_n=%g so_p=%g pot_n=%g pot_p=%g mag=%g (rel=%g)\n", *((double*)(i->data)), ((avg_n+pot_n)-(avg_p+pot_p)), ((so_n+pot_n)-(so_p+pot_p)), (avg_n-relative_to), (avg_p-relative_to), (so_n-relative_to), (so_p-relative_to), (pot_n), (pot_p), complex_magnitude(gsl_spline_eval(re, *((double*)(i->data)), re_acc), gsl_spline_eval(im, *((double*)(i->data)), im_acc)), relative_to);
    alpha += complex_magnitude(gsl_spline_eval(re, *((double*)(i->data)), re_acc), gsl_spline_eval(im, *((double*)(i->data)), im_acc))
      *(
  ((so_n+pot_n)-(so_p+pot_p))*(1.0/((so_n)*(so_n))+1.0/((so_p)*(so_p)))
  -((avg_n+pot_n)-(avg_p+pot_p))*(1.0/((avg_n)*(avg_n))+1.0/((avg_p)*(avg_p)))
  );
  }
  /*De-allocate structures.*/
  /*HH*/
  gsl_spline_free(hh);
  gsl_interp_accel_free(hh_acc);
  /*LH*/
  gsl_spline_free(lh);
  gsl_interp_accel_free(lh_acc);
  /*SO*/
  gsl_spline_free(so);
  gsl_interp_accel_free(so_acc);
  /*CB*/
  gsl_spline_free(cb);
  gsl_interp_accel_free(cb_acc);
  /*pot*/
  gsl_spline_free(pot);
  gsl_interp_accel_free(pot_acc);
  /*wavefunc*/
  gsl_spline_free(re);
  gsl_interp_accel_free(re_acc);
  /*HH*/
  gsl_spline_free(im);
  gsl_interp_accel_free(im_acc);
  free(repacked_x);
  free(repacked_re);
  free(repacked_im);
  free(avg);
  /*NOTE that the one-half from the summation has been absorbed into here.*/
  /*NOTE that the minus sign comes from the definitions of
    delta_gamma_x and S_n.  See notebook on April 28, 2009*/
  /*Don't forget to convert nm->m (from the wavefunction)*/
  //alpha *= hbar*hbar*22/12/m_e_0*1e9;  /*22eV, from the paper*/
  alpha *= hbar*hbar*Ep/12/m_e_0*1e9;  /*22eV, from the paper*/
  return alpha;
}

struct integration_param {
  gsl_spline *re;
  gsl_interp_accel *re_accel;
  gsl_spline *im;
  gsl_interp_accel *im_accel;

  gsl_spline *hh;
  gsl_interp_accel *hh_accel;
  gsl_spline *lh;
  gsl_interp_accel *lh_accel;
  gsl_spline *so;
  gsl_interp_accel *so_accel;
  gsl_spline *cb;
  gsl_interp_accel *cb_accel;

  gsl_spline *avg;
  gsl_interp_accel *avg_accel;
  
  gsl_spline *pot;
  gsl_interp_accel *pot_accel;
};

struct integration_param_2 {
  gsl_spline *re;
  gsl_interp_accel *re_accel;
  gsl_spline *im;
  gsl_interp_accel *im_accel;

  gsl_spline *hh;
  gsl_interp_accel *hh_accel;
  gsl_spline *lh;
  gsl_interp_accel *lh_accel;
  gsl_spline *so;
  gsl_interp_accel *so_accel;
  gsl_spline *cb;
  gsl_interp_accel *cb_accel;

  gsl_spline *avg;
  gsl_interp_accel *avg_accel;
  
  gsl_spline *pot;
  gsl_interp_accel *pot_accel;

  double Ep;
  gsl_spline *Ep_z;
  gsl_interp_accel *Ep_z_accel;
};

__always_inline double band_terms(double x, void *params) {
  //register double pot = gsl_spline_eval(((struct integration_param*)params)->pot, x, ((struct integration_param*)params)->pot_accel);
  double avg = (gsl_spline_eval(((struct integration_param*)params)->avg, x, ((struct integration_param*)params)->avg_accel));
  double so = (gsl_spline_eval(((struct integration_param*)params)->so, x, ((struct integration_param*)params)->so_accel));
  return (1.0/(so*so) - 1.0/(avg*avg));
}

/*E_F is at 0*/
double integrand(double x, void *params) {
  register double re = gsl_spline_eval(((struct integration_param*)params)->re, x, ((struct integration_param*)params)->re_accel);
  register double im = gsl_spline_eval(((struct integration_param*)params)->im, x, ((struct integration_param*)params)->im_accel);
  register double ef = gsl_spline_eval_deriv(((struct integration_param*)params)->pot, x, ((struct integration_param*)params)->pot_accel);
  return ef*(re*re + im*im)*band_terms(x, params);
}

__always_inline double band_terms_del_psi2_method(double x, void *params) {
  //register double pot = gsl_spline_eval(((struct integration_param_2*)params)->pot, x, ((struct integration_param_2*)params)->pot_accel);
  double avg = (gsl_spline_eval(((struct integration_param_2*)params)->avg, x, ((struct integration_param_2*)params)->avg_accel));
  double so = (gsl_spline_eval(((struct integration_param_2*)params)->so, x, ((struct integration_param_2*)params)->so_accel));
  return (1.0/so - 1.0/avg);
}

/*d/dz(re^2 + im^2) = 2*(re*d/dz re + im*d/dz im)*/
__always_inline double del_psi2(double x, gsl_spline *re, gsl_interp_accel *re_accel, gsl_spline *im, gsl_interp_accel *im_accel) {
  //fprintf(stderr, "%g: %g->%g %g->%g\n", x, gsl_spline_eval(re, x, re_accel), gsl_spline_eval_deriv(re, x, re_accel), gsl_spline_eval(im, x, im_accel), gsl_spline_eval_deriv(im, x, im_accel));
  return 2*(gsl_spline_eval(re, x, re_accel)*gsl_spline_eval_deriv(re, x, re_accel) + gsl_spline_eval(im, x, im_accel)*gsl_spline_eval_deriv(im, x, im_accel));
}

/*E_F is at 0*/
double integrand_del_psi2_method(double x, void *params) {
  if(((struct integration_param_2*)params)->Ep_z == NULL) {
    return (((struct integration_param_2*)params)->Ep) 
      *del_psi2(x, ((struct integration_param_2*)params)->re, ((struct integration_param_2*)params)->re_accel, ((struct integration_param_2*)params)->im, ((struct integration_param_2*)params)->im_accel)
      *band_terms_del_psi2_method(x, params);
  }else{
    return gsl_spline_eval(((struct integration_param_2*)params)->Ep_z, x, ((struct integration_param_2*)params)->Ep_z_accel)
      *del_psi2(x, ((struct integration_param_2*)params)->re, ((struct integration_param_2*)params)->re_accel, ((struct integration_param_2*)params)->im, ((struct integration_param_2*)params)->im_accel)
      *band_terms_del_psi2_method(x, params);
  }
}

double alpha_electric_term_mstar_nostrain(const struct mstar_wavefunction_1d *wf, const struct band *hh_band, const struct band *lh_band, const struct band *so_band, const struct band *cb_band, const struct potential *potential, struct alpha_integrand_verificator *vfp, double *integration_error, double abserr, double relerr, gboolean no_integrate, double Ep) {
  /*HH*/
  gsl_interp_accel *hh_acc = gsl_interp_accel_alloc();
  gsl_spline *hh = gsl_spline_alloc(gsl_interp_akima, hh_band->N);
  gsl_spline_init(hh, hh_band->x, hh_band->bandedge, hh_band->N);
  /*LH*/
  gsl_interp_accel *lh_acc = gsl_interp_accel_alloc();
  gsl_spline *lh = gsl_spline_alloc(gsl_interp_akima, lh_band->N);
  gsl_spline_init(lh, lh_band->x, lh_band->bandedge, lh_band->N);
  /*SO*/
  gsl_interp_accel *so_acc = gsl_interp_accel_alloc();
  gsl_spline *so = gsl_spline_alloc(gsl_interp_akima, so_band->N);
  gsl_spline_init(so, so_band->x, so_band->bandedge, so_band->N);
  /*CB*/
  gsl_interp_accel *cb_acc = gsl_interp_accel_alloc();
  gsl_spline *cb = gsl_spline_alloc(gsl_interp_akima, cb_band->N);
  gsl_spline_init(cb, cb_band->x, cb_band->bandedge, cb_band->N);
  /*Potential*/
  gsl_interp_accel *pot_acc = gsl_interp_accel_alloc();
  gsl_spline *pot = gsl_spline_alloc(gsl_interp_akima, potential->N);
  gsl_spline_init(pot, potential->x, potential->potential, potential->N);
  /*Get the average of the heavy and light holes for now*/
  double *arrays[] = {hh_band->bandedge, lh_band->bandedge};
  double *avg_array = average_arrays(arrays, 2, hh_band->N);
  gsl_interp_accel *avg_acc = gsl_interp_accel_alloc();
  gsl_spline *avg = gsl_spline_alloc(gsl_interp_akima, hh_band->N);
  gsl_spline_init(avg, hh_band->x, avg_array, hh_band->N);
  /*Wavefunction*/
  double* repacked_x;
  double* repacked_re;
  double* repacked_im;
  //for(unsigned i=0; i<wf->N; i++) fprintf(stderr, "wf.x=%g wf.re=%g\n", wf->points[i].x, wf->points[i].re);
  mstar_wavefunction_1d_repack(wf, &repacked_x, &repacked_re, &repacked_im, 1);
  //for(unsigned i=0; i<wf->N; i++) fprintf(stderr, "->repacked_x=%g repacked_re=%g\n", repacked_x[i], repacked_re[i]);
  gsl_interp_accel *re_acc = gsl_interp_accel_alloc();
  gsl_spline *re = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(re, repacked_x, repacked_re, wf->N);
  gsl_interp_accel *im_acc = gsl_interp_accel_alloc();
  gsl_spline *im = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(im, repacked_x, repacked_im, wf->N);
  /*DO THE INTEGRATION*/
  struct integration_param parms = {
    .re = re,
    .re_accel = re_acc,
    .im = im,
    .im_accel = im_acc,
    .hh = hh,
    .hh_accel = hh_acc,
    .lh = lh,
    .lh_accel = lh_acc,
    .so = so,
    .so_accel = so_acc,
    .cb = cb,
    .cb_accel = cb_acc,
    .avg = avg,
    .avg_accel = avg_acc,
    .pot = pot,
    .pot_accel = pot_acc,
  };
  gsl_function gsl_func = {
    .function = &integrand,
    .params = &parms
  };
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(MAX_GSL_INTEGRATION_INTERVALS);
  double alpha;
  int status=0;
  if(!no_integrate) status = gsl_integration_qag(&gsl_func, repacked_x[0], repacked_x[wf->N-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w, &alpha, integration_error);
  if(status != 0) {
    fprintf(stderr, "alpha_electric_field_term_mstar_nostrain: ERROR occurred while integrating: %d (%s); returned %g w/error %g.  Retrying with GAUSS31\n", status, gsl_strerror(status), alpha, *integration_error);
    status = gsl_integration_qag(&gsl_func, repacked_x[0], repacked_x[wf->N-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS31, w, &alpha, integration_error);
    if(status != 0) {
      fprintf(stderr, "\tFAILED: %d (%s); returned %g w/error %g   Retrying with GAUSS15\n", status, gsl_strerror(status), alpha, *integration_error);
      status = gsl_integration_qag(&gsl_func, repacked_x[0], repacked_x[wf->N-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS15, w, &alpha, integration_error);
      if(status != 0) {
  size_t neval;
  fprintf(stderr, "\t\tFAILEd: %d (%s); returned %g w/error %g  Retrying with non-adaptive Gauss-Kronrod\n", status, gsl_strerror(status), alpha, *integration_error);
  status = gsl_integration_qng (&gsl_func, repacked_x[0], repacked_x[wf->N-1], abserr, relerr, &alpha, integration_error, &neval);
  if(status != 0) {
    fprintf(stderr, "\t\t\tFAILED : %d (%s); returned %g w/error %g (%ld evals)\n", status, gsl_strerror(status), alpha, *integration_error, (long int)neval);
  }else{
    fprintf(stderr, "\t\t\tSUCCESS: got %g w/error %g (%ld evals)\n", alpha, *integration_error, (long int)neval);
  }
      }else{
  fprintf(stderr, "\t\tSUCCESS: got %g w/error %g\n", alpha, *integration_error);
      }
    }else{
      fprintf(stderr, "\tSUCCESS: got %g w/error %g\n", alpha, *integration_error);
    }
  }
  gsl_integration_workspace_free(w);
  /*Save off the verification info, if needed.*/
  if(vfp != NULL) {
    vfp->N = 5000;
    vfp->x = malloc((vfp->N)*sizeof(double));
    vfp->prob_density = malloc((vfp->N)*sizeof(double));
    vfp->integrand = malloc((vfp->N)*sizeof(double));
    vfp->pot_deriv = malloc((vfp->N)*sizeof(double));
    vfp->pot = malloc((vfp->N)*sizeof(double));
    vfp->iface_equiv = malloc((vfp->N)*sizeof(double));
    vfp->sob = malloc((vfp->N)*sizeof(double));
    vfp->avg = malloc((vfp->N)*sizeof(double));
    vfp->oo_sob = malloc((vfp->N)*sizeof(double));
    vfp->oo_avg = malloc((vfp->N)*sizeof(double));
    for(int i=0; i<(vfp->N); i++) {
      (vfp->x)[i] = (potential->x)[0] + ((potential->x)[potential->N-1] - (potential->x)[0])/vfp->N*i;
      (vfp->integrand)[i] = integrand((vfp->x)[i], &parms);
      (vfp->prob_density)[i] = complex_magnitude(gsl_spline_eval(re, (vfp->x)[i], re_acc), gsl_spline_eval(im, (vfp->x)[i], im_acc));
      (vfp->pot_deriv)[i] = gsl_spline_eval_deriv(pot, (vfp->x)[i], pot_acc);
      (vfp->pot)[i] = gsl_spline_eval(pot, (vfp->x)[i], pot_acc);
      (vfp->iface_equiv)[i] = band_terms((vfp->x)[i], &parms);
      (vfp->sob)[i] = gsl_spline_eval(so, (vfp->x)[i], so_acc);
      (vfp->avg)[i] = gsl_spline_eval(avg, (vfp->x)[i], avg_acc);
      (vfp->oo_sob)[i] = 1.0/(vfp->sob)[i];
      (vfp->oo_avg)[i] = 1.0/(vfp->avg)[i];
    }
  }
  /*De-allocate structures.*/
  /*HH*/
  gsl_spline_free(hh);
  gsl_interp_accel_free(hh_acc);
  /*LH*/
  gsl_spline_free(lh);
  gsl_interp_accel_free(lh_acc);
  /*SO*/
  gsl_spline_free(so);
  gsl_interp_accel_free(so_acc);
  /*CB*/
  gsl_spline_free(cb);
  gsl_interp_accel_free(cb_acc);
  /*pot*/
  gsl_spline_free(pot);
  gsl_interp_accel_free(pot_acc);
  /*wavefunc*/
  gsl_spline_free(re);
  gsl_interp_accel_free(re_acc);
  /*HH*/
  gsl_spline_free(im);
  gsl_interp_accel_free(im_acc);
  free(repacked_x);
  free(repacked_re);
  free(repacked_im);
  free(avg);
  /*The factor of - is required to convert the volts to eV*/
  /*Don't forget to convert nm->m (from the wavefunction)*/
  *integration_error *= -hbar*hbar*Ep/6*oo_m_e_0*1e9;
  //alpha *= hbar*hbar*22/6*oo_m_e_0*1e9;
  alpha *= hbar*hbar*Ep/6*oo_m_e_0*1e9;
  return alpha;
}


double alpha_from_del_psi2_nostrain(const struct mstar_wavefunction_1d *wf, const struct band *hh_band, const struct band *lh_band, const struct band *so_band, const struct band *cb_band, const struct potential *potential, struct alpha_integrand_verificator_2 *vfp2, double *integration_error, double abserr, double relerr, gboolean no_integrate, double Ep) {
  /*HH*/
  gsl_interp_accel *hh_acc = gsl_interp_accel_alloc();
  gsl_spline *hh = gsl_spline_alloc(gsl_interp_akima, hh_band->N);
  gsl_spline_init(hh, hh_band->x, hh_band->bandedge, hh_band->N);
  /*LH*/
  gsl_interp_accel *lh_acc = gsl_interp_accel_alloc();
  gsl_spline *lh = gsl_spline_alloc(gsl_interp_akima, lh_band->N);
  gsl_spline_init(lh, lh_band->x, lh_band->bandedge, lh_band->N);
  /*SO*/
  gsl_interp_accel *so_acc = gsl_interp_accel_alloc();
  gsl_spline *so = gsl_spline_alloc(gsl_interp_akima, so_band->N);
  gsl_spline_init(so, so_band->x, so_band->bandedge, so_band->N);
  /*CB*/
  gsl_interp_accel *cb_acc = gsl_interp_accel_alloc();
  gsl_spline *cb = gsl_spline_alloc(gsl_interp_akima, cb_band->N);
  gsl_spline_init(cb, cb_band->x, cb_band->bandedge, cb_band->N);
  /*Potential*/
  gsl_interp_accel *pot_acc = gsl_interp_accel_alloc();
  gsl_spline *pot = gsl_spline_alloc(gsl_interp_akima, potential->N);
  gsl_spline_init(pot, potential->x, potential->potential, potential->N);
  /*Get the average of the heavy and light holes for now*/
  double *arrays[] = {hh_band->bandedge, lh_band->bandedge};
  double *avg_array = average_arrays(arrays, 2, hh_band->N);
  gsl_interp_accel *avg_acc = gsl_interp_accel_alloc();
  gsl_spline *avg = gsl_spline_alloc(gsl_interp_akima, hh_band->N);
  gsl_spline_init(avg, hh_band->x, avg_array, hh_band->N);
  /*Wavefunction*/
  double* repacked_x;
  double* repacked_re;
  double* repacked_im;
  mstar_wavefunction_1d_repack(wf, &repacked_x, &repacked_re, &repacked_im, 1);
  gsl_interp_accel *re_acc = gsl_interp_accel_alloc();
  /*gsl_spline *re = gsl_spline_alloc(gsl_interp_cspline, wf->N);*/
  gsl_spline *re = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(re, repacked_x, repacked_re, wf->N);
  gsl_interp_accel *im_acc = gsl_interp_accel_alloc();
  /*gsl_spline *im = gsl_spline_alloc(gsl_interp_cspline, wf->N);*/
  gsl_spline *im = gsl_spline_alloc(gsl_interp_akima, wf->N);
  gsl_spline_init(im, repacked_x, repacked_im, wf->N);
  /*DO THE INTEGRATION*/
  struct integration_param_2 parms = {
    .re = re,
    .re_accel = re_acc,
    .im = im,
    .im_accel = im_acc,
    .hh = hh,
    .hh_accel = hh_acc,
    .lh = lh,
    .lh_accel = lh_acc,
    .so = so,
    .so_accel = so_acc,
    .cb = cb,
    .cb_accel = cb_acc,
    .avg = avg,
    .avg_accel = avg_acc,
    .pot = pot,
    .pot_accel = pot_acc,
    .Ep = Ep,
    .Ep_z = NULL,
    .Ep_z_accel = NULL
  };
  gsl_function gsl_func = {
    .function = & integrand_del_psi2_method,
    .params = &parms
  };
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(MAX_GSL_INTEGRATION_INTERVALS);
  double alpha;
  int status=0;
  if(!no_integrate) status = gsl_integration_qag(&gsl_func, repacked_x[0], repacked_x[wf->N-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w, &alpha, integration_error);
  if(status != 0) {
    fprintf(stderr, "alpha_calc::alpha_from_del_psi2_nostrain: ERROR occurred while integrating: %d (%s)\n", status, gsl_strerror(status));
  }
  gsl_integration_workspace_free(w);
  /*Save off the verification info, if needed.*/
  if(vfp2 != NULL) {
    vfp2->N = 5000;
    vfp2->x = malloc((vfp2->N)*sizeof(double));
    vfp2->prob_density = malloc((vfp2->N)*sizeof(double));
    vfp2->d_prob_density = malloc((vfp2->N)*sizeof(double));
    vfp2->integrand = malloc((vfp2->N)*sizeof(double));
    vfp2->sob = malloc((vfp2->N)*sizeof(double));
    vfp2->avg = malloc((vfp2->N)*sizeof(double));
    vfp2->oo_sob = malloc((vfp2->N)*sizeof(double));
    vfp2->oo_avg = malloc((vfp2->N)*sizeof(double));
    vfp2->band_term = malloc((vfp2->N)*sizeof(double));
    for(int i=0; i<(vfp2->N); i++) {
      (vfp2->x)[i] = (potential->x)[0] + ((potential->x)[potential->N-1] - (potential->x)[0])/vfp2->N*i;
      (vfp2->integrand)[i] = integrand_del_psi2_method((vfp2->x)[i], &parms);
      (vfp2->prob_density)[i] = complex_magnitude(gsl_spline_eval(re, (vfp2->x)[i], re_acc), gsl_spline_eval(im, (vfp2->x)[i], im_acc));
      (vfp2->d_prob_density)[i] = del_psi2((vfp2->x)[i], re, re_acc, im, im_acc);
      (vfp2->sob)[i] = gsl_spline_eval(so, (vfp2->x)[i], so_acc);
      (vfp2->avg)[i] = gsl_spline_eval(avg, (vfp2->x)[i], avg_acc);
      (vfp2->oo_sob)[i] = 1.0/(vfp2->sob)[i];
      (vfp2->oo_avg)[i] = 1.0/(vfp2->avg)[i];
      (vfp2->band_term)[i] = band_terms_del_psi2_method((vfp2->x)[i], &parms);
    }
  }
  /*De-allocate structures.*/
  /*HH*/
  gsl_spline_free(hh); 
  gsl_interp_accel_free(hh_acc);
  /*LH*/
  gsl_spline_free(lh);
  gsl_interp_accel_free(lh_acc);
  /*SO*/
  gsl_spline_free(so);
  gsl_interp_accel_free(so_acc);
  /*CB*/
  gsl_spline_free(cb);
  gsl_interp_accel_free(cb_acc);
  /*pot*/
  gsl_spline_free(pot);
  gsl_interp_accel_free(pot_acc);
  /*wavefunc*/
  gsl_spline_free(re);
  gsl_interp_accel_free(re_acc);
  /*HH*/
  gsl_spline_free(im);
  gsl_interp_accel_free(im_acc);
  free(repacked_x);
  free(repacked_re);
  free(repacked_im);
  /*The factor of - is required to convert the volts to eV*/
  /*Don't forget to convert nm->m (from the wavefunction)*/
  *integration_error *= -hbar*hbar*6*oo_m_e_0*1e9;
  //alpha *= hbar*hbar*22/6*oo_m_e_0*1e9;
  alpha *= hbar*hbar*6*oo_m_e_0*1e9;
  return alpha;
}

GError* alpha_write_verificator_data(GFile* f, const struct alpha_integrand_verificator* vf) {
  GError *gerror = NULL;
  GFileOutputStream *os = g_file_replace(f, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerror);
  if(os == NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, gerror, "ERROR when opening alpha_efield instrumentation file for writing");
    return e;
  }
  GString *str = g_string_new("");
  GDataOutputStream *gdos = g_data_output_stream_new((GOutputStream*)os);
  g_string_printf(str, "#x(nm)\tprob_density(1/nm)\tintegrand(1/(ev^2 nm))\tpot_deriv(v/nm)\tband_terms(1/ev^2)\tsplit-off (eV)\taverged HH,LH (eV)\t1/split-off (1/eV)\t1/averaged (1/eV)\tPotential (eV)\n");
  g_data_output_stream_put_string(gdos, str->str, NULL, &gerror);
  if(gerror != NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, gerror, "ERROR when writing alpha_efield instrumentation file header");
    g_object_unref(os);
    g_object_unref(gdos);
    g_string_free(str, TRUE);
    return e;
  }
  int N;
  N = vf->N;
  for(int i=0; i<N; i++) {
    g_string_printf(str, "%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", (vf->x)[i], (vf->prob_density)[i], (vf->integrand)[i], (vf->pot_deriv)[i], (vf->iface_equiv)[i], (vf->sob)[i], (vf->avg)[i], (vf->oo_sob)[i], (vf->oo_avg)[i], (vf->pot)[i]);
    g_data_output_stream_put_string(gdos, str->str, NULL, &gerror);
    if(gerror != NULL) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, gerror, "ERROR when writing alpha_efield instrumentation file header (Missing line would have been:%s)", str->str);
      g_object_unref(os);
      g_object_unref(gdos);
      g_string_free(str, TRUE);
      return e;
    }
  }
  g_output_stream_close(G_OUTPUT_STREAM(gdos), NULL, &gerror);
  if(gerror != NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, gerror, "ERROR when closing alpha_efield instrumentation output stream.\n");
  }
  g_object_unref(gdos);
  g_object_unref(os);
  g_string_free(str, TRUE);
  return NULL;
}
GError* alpha_write_verificator2_data(GFile* f, const struct alpha_integrand_verificator_2* vf2) {
  GError *gerror = NULL;
  GFileOutputStream *os = g_file_replace(f, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerror);
  if(os == NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, gerror, "ERROR when opening method 2 instrumentation file for writing");
    return e;
  }
  GString *str = g_string_new("");
  GDataOutputStream *gdos = g_data_output_stream_new((GOutputStream*)os);
  g_string_printf(str, "#x(nm)\tprob_density(1/nm)\td_dz prob_density\tintegrand(1/(ev^2 nm))\tsplit-off (eV)\taverged HH,LH (eV)\t1/split-off (1/eV)\t1/averaged (1/eV)\tband terms (1eV)\n");
  g_data_output_stream_put_string(gdos, str->str, NULL, &gerror);
  if(gerror != NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, gerror, "ERROR when writing method 2 instrumentation file header");
    g_object_unref(os);
    g_object_unref(gdos);
    g_string_free(str, TRUE);
    return e;
  }
  int N;
  N = vf2->N;
  for(int i=0; i<N; i++) {
    g_string_printf(str, "%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", (vf2->x)[i], (vf2->prob_density)[i], (vf2->d_prob_density)[i], (vf2->integrand)[i], (vf2->sob)[i], (vf2->avg)[i], (vf2->oo_sob)[i], (vf2->oo_avg)[i], (vf2->band_term)[i]);
    g_data_output_stream_put_string(gdos, str->str, NULL, &gerror);
    if(gerror != NULL) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, gerror, "ERROR when writing method 2 instrumentation file header (Missing line would have been:%s)", str->str);
      g_object_unref(os);
      g_object_unref(gdos);
      g_string_free(str, TRUE);
      return e;
    }
  }
  g_output_stream_close(G_OUTPUT_STREAM(gdos), NULL, &gerror);
  if(gerror != NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, gerror, "ERROR when closing method 2 instrumentation output stream.\n");
  }
  g_object_unref(gdos);
  g_object_unref(os);
  g_string_free(str, TRUE);
  return NULL;
}

void alpha_free_verificator_data(struct alpha_integrand_verificator *vf) {
  free(vf->x);
  free(vf->prob_density);
  free(vf->integrand);
  free(vf->pot_deriv);
  free(vf->iface_equiv);
  free(vf->sob);
  free(vf->avg);
  free(vf->oo_sob);
  free(vf->oo_avg);
  free(vf->pot);
}

void alpha_free_verificator2_data(struct alpha_integrand_verificator_2* vf2) {
  free(vf2->x);
  free(vf2->prob_density);
  free(vf2->d_prob_density);
  free(vf2->band_term);
  free(vf2->sob);
  free(vf2->avg);
  free(vf2->oo_sob);
  free(vf2->oo_avg);
  free(vf2->integrand);
}