#include <iostream>

const float cVac=2.998e-1;
const float km2m=1e3;

float ntimexpASlib(float z1, float z2, float dtot){
  z1 /= 1000;  z2 /= 1000;  dtot /= 1000;  float A=1.78;  float B=A*-0.24;
  const float C=16.; //was 16. but must be shifted since units are km
  float zmax=z2;
  float zmin=z1;
  if(zmin>zmax){
    zmin=z2;
    zmax=z1;
  }
  float ctheta;
  if(dtot>0)
    ctheta=fabs(z2-z1)/dtot;
  else
    return 0;
  float ftimexp=dtot/cVac; //was =0
  if(ctheta!=0){
    if(zmax<0. && zmin<0.)
      ftimexp=(A*(zmax-zmin)+(B/C)*(exp(C*zmax)-exp(C*zmin)))/(cVac * ctheta);
    if(zmax>=0. && zmin<0.)
      ftimexp=(A*(0.-zmin)+(B/C)*(exp(C*0.)-exp(C*zmin))+zmax)/(cVac * ctheta);
    if(zmax>=0. && zmin>=0.)
      ftimexp=(zmax-zmin)/(cVac * ctheta);
    if(zmax<0. && zmin>0.)
      cout<<" CONDITION NOT COVERED IN NTIMEZ!!!"<<endl;
  }
  return ftimexp*km2m; //km->meters
}

double InvTimeExp(double time, double z1, double z2) {

  z1 /= km2m;
  z2 /= km2m;
  time /= km2m;

  double zmax = z1;
  double zmin = z2;

  if (zmin>zmax) {

    zmin = z1;
    zmax = z2;

  }

  double A = 1.78;
  double B = A * (-.24);
  double C = 16;

  double dist = (time * cVac * fabs(zmax-zmin) ) /
    ( A * (zmax-zmin) + (B/C)* ( exp(C*zmax) - exp(C*zmin)) );

  return (dist*km2m);

}

double SQ(double a) {

  return(a*a);

}

double Dist3d(double x1, double y1, double z1,
	      double x2, double y2, double z2) {

  return ( sqrt( SQ(x2-x1) + SQ(y2-y1) + SQ(z2-z1) ) );

}