/* deMunck's multilayer model. */
/* 
   See J.C. Demunck and Maria J. Peters,
   "A Fast Method to Compute the Potential 
   in the Multisphere Model", IEEE Trans. 
   Biomed. Eng. , Vo 40. N. 11, November 1993.
   pp 1168-1174
   All Units are MKS, dammit.
 */
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include <limits.h>
/*#include <values.h>*/
#include "matrixmalloc.h"
#include <unistd.h>
#include <getopt.h>
#define FOUR_PI 12.5664f
#define Pi 3.14159265358979
#define pi Pi
#include <signal.h>


int I,J,K;
int start_i=1,start_j=1,start_k=1;
int curr_i,curr_j,curr_k;

int minimum_terms=50;

int slow=0,singleslice=0,sensorlist=0;
int isotropic =1; 

double rad,mono=0.0f;
int filemode; 
double *r, *eps, *eta;
int layers;
double Mrad =0.0f,  Mtan =0.0f ;
double  dx=1.0f,dy=1.0f,dz=1.0f,rbottom=0.0f,rbottom2=0.0f;
double sx=0.0f,sy=0.0f,sz=0.0f;
double tolerance =0.001f; 
double *Rbottom, *Rright; 
double * legendre;
int rbsize; 
void oopsie(int signum) { 
  FILE * disasterfile; 
  fprintf(stderr, "Caught signal %d, while on voxel %d %d %d.",
	 signum,curr_i,curr_j,curr_k); 
  disasterfile=fopen("demunck.disaster","w");
  fprintf(disasterfile, "Caught signal %d, while on voxel %d %d %d.",
	 signum,curr_i,curr_j,curr_k); 
  fclose(disasterfile); 
  exit(-22); 
} 

#if 0
int ijkton(int i, int j, int k){
  return i+I*(j-1)+I*J*(k-1);
}
#endif

double nu(int j,int n)
{ /* equation 14 */
  return
    isotropic?
    n:
  ((-1+sqrt(1+
	    4*n*(n+1)*eta[j]/eps[j]))
   /2);
}

void smatrix(int sign, int j, double r_a, double r_b,int n, double out[3][3]) {
  /* equations 30 and 29 */
  double nu_val,nu2;

  nu_val=nu(j,n);
  nu2= 2*nu_val+1;
  if (sign == -1) {
    out[1][1] = nu_val*r_a/(nu2*r_b);
    out[1][2] = -r_a/(nu2*eps[j]);
    out[2][1] = -(nu_val+1)*nu_val*eps[j]/(nu2*r_b);
    out[2][2] = (nu_val+1)/nu2;
  }  
  if (sign == 1) {
    out[1][1] = (nu_val+1)/nu2;
    out[1][2] = r_b/(nu2*eps[j]);
    out[2][1] = (nu_val+1)*nu_val*eps[j]/(nu2*r_a);
    out[2][2] = nu_val*r_b/(nu2*r_a);
  }

}

double beta(int j) {
  return sqrt(eta[j]*eps[j]);
}

double alpha(int j ) {
  return sqrt(eta[j]/eps[j]);
}

void matrix_prod_2x2(double a[3][3], double b[3][3], double prod[3][3]){

  prod[1][1] = a[1][1]*b[1][1]+a[1][2]*b[2][1];
  prod[1][2] = a[1][1]*b[1][2]+a[1][2]*b[2][2];
  prod[2][1] = a[2][1]*b[1][1]+a[2][2]*b[2][1];
  prod[2][2] = a[2][1]*b[1][2]+a[2][2]*b[2][2];

}

void matrix_prod(double** a, double** b, double** prod){

  prod[1][1] = a[1][1]*b[1][1]+a[1][2]*b[2][1];
  prod[1][2] = a[1][1]*b[1][2]+a[1][2]*b[2][2];
  prod[2][1] = a[2][1]*b[1][1]+a[2][2]*b[2][1];
  prod[2][2] = a[2][1]*b[1][2]+a[2][2]*b[2][2];

}

void clear_2x2 (double **matrix){

  free(matrix[1]);

  free(matrix);
}

void matrix_sum_2x2(double a[3][3], double b[3][3]){
  a[1][1]+=b[1][1];
  a[1][2]+=b[1][2];
  a[2][2]+=b[2][2];
  a[2][1]+=b[2][1];
  
}

void copy_2x2 (double a[3][3], double b[3][3]){
  a[1][1]=b[1][1];
  a[1][2]=b[1][2];
  a[2][1]=b[2][1];
  a[2][2]=b[2][2];
}

void matrix_scalar_2x2(double a[3][3], double b){
  a[1][1]*=b;
  a[1][2]*=b;
  a[2][1]*=b;
  a[2][2]*=b;
}
double lambda(int j,double nu_){/* equation 33 */
  return pow(r[j+1]/r[j],nu_);
}

void u_matrix(int n,int j, double r_a, double r_b, double out[3][3]){
  /* equation 31 */
  double sminus[3][3];
  double splus[3][3];
  double tempscalar,nu_;

  if (r_b>0.0f) { 
    
    tempscalar=pow(r_a/r_b,-2*nu(j,n)-2);
    
    smatrix(1,j,r_a,r_b,n,splus);/* r_a r_b*/
    smatrix(-1,j,r_a,r_b,n,sminus);
    matrix_scalar_2x2(sminus,tempscalar);
  
    matrix_sum_2x2(splus,sminus);
    copy_2x2(out,splus);
 
  }else {
    /* the problem for the deep cases lies HERE. */

    nu_=nu(layers,n);
    out[1][1]=0.0f; 
    out[2][1]=0.0f;
    out[1][2]=1/(eps[layers]*(1+2*nu_)); 
    out[2][2]=nu_/(r_a*(1+2*nu_)); 
    if (!n) {
      out[1][2]=-1/(r_a*eps[layers]);
      out[2][2]=1/(r_a*r_a);
    }
  }

}


double R(int n,double r_o, double r_e,
	 int Jo, int Je) {
  /* equation 32 */
  double nu_, nu_lJo, nu_lJe;
  int j,lJe,lJo;
  double lr_o,lr_e;
  double Rout;
  double mat_[3][3]; 
  double mat_a[3][3]; 
  double mat_b[3][3];
  double mat_c[3][3];                    
  double mat_d[3][3];
  /* suspect situations : 
     when there is a border between source and point. 
     R() is symmetric when there aren't any borders. */

  if (r_o <= r_e) { 
    lr_o=r_o;
    lr_e=r_e;
    lJe=Je;
    lJo=Jo;
  } else {
    /* modify things for electrode deeper than dipole */
    lr_o=r_e;
    lr_e=r_o;
    lJe=Jo;
    lJo=Je;
  }
  nu_lJo=nu(lJo,n);
  nu_lJe=nu(lJe,n);

  if (lJe > 1) { 
    u_matrix(n,1,r[1],r[2],mat_c);
    for (j=2; j<=lJe-1; j++) {
      
      u_matrix(n,j,r[j],r[j+1],mat_);

      copy_2x2(mat_d,mat_c);
      matrix_prod_2x2(mat_d,mat_,mat_c);
    }
    u_matrix(n,lJe,r[lJe],lr_e,mat_a);

    matrix_prod_2x2(mat_c,mat_a,mat_d);
    Rout= mat_d[2][2];
  } else { 
    if (!lJe) 
      Rout=1.0f;
    else { /* lje==1 */
      u_matrix(n,1,r[1],lr_e,mat_a);
      Rout= mat_a[2][2];
    }
  }

  u_matrix(n,lJo,lr_o,r[lJo+1],mat_c);
    
    for (j=lJo+1; j<=layers; j++) {
      copy_2x2(mat_d,mat_c);
      u_matrix(n,j,r[j],r[j+1],mat_);
      matrix_prod_2x2(mat_d,mat_,mat_c);
    }
    
    Rout*=
      mat_c[1][2];

    if (n > rbsize) {
      printf("-"); fflush(NULL);
      rbsize +=200;
      Rbottom = (double*)realloc(Rbottom,sizeof(double)*(rbsize+1));
      legendre = (double*)realloc(legendre,sizeof(double)*(rbsize+2));
      for (j=rbsize-200; j <= rbsize; j++)
	Rbottom[j]=0.0f; 
      
    }
    if (Rbottom[n] == 0.0f) {  
      u_matrix(n,1,r[1],r[2],mat_c);
      for (j=2; j<=layers; j++) {
	copy_2x2(mat_d,mat_c);
	u_matrix(n,j,r[j],r[j+1],mat_);
	matrix_prod_2x2(mat_d,mat_,mat_c);
      }
      Rbottom[n] = -mat_c[2][2];
    }
    /* this ought to be very suspect: */
    Rout/=lr_e*lr_e*Rbottom[n];
  //  Rout/=r[1]*r[1]*Rbottom[n];

  Rout *=
    pow(lr_o/r[lJo],nu_lJo)/
    pow(lr_e/r[lJe],nu_lJe);
  for (j=lJe; j<lJo; j++) {
    nu_=nu(j,n);
    Rout *=lambda(j,nu_);
    
  }

  return Rout;

}

double Rprime(int n,double r_o, double r_e,
	 int Jo, int Je) {
  /* equation 60 */
  int j,lJe,lJo,deep;
  double nu_, nu_lJo, nu_lJe;
  double lr_o,lr_e;
  double Rpr;
  double mat_[3][3]; 
  double mat_a[3][3]; 
  double mat_b[3][3];
  double mat_c[3][3];   
  double mat_d[3][3];


  if (r_o <= r_e) { 
    deep =0;
    lr_o=r_o;
    lr_e=r_e;
    lJe=Je;
    lJo=Jo;
  } else {
    /* modify things for electrode deeper than dipole */
    deep=1;
    lr_o=r_e;
    lr_e=r_o;
    lJe=Jo;
    lJo=Je;
  }

  nu_lJo=nu(lJo,n);
  nu_lJe=nu(lJe,n);
  /* outermost shell down to first of r_e and r_o */
  /* can be automated only for the deep case. */
  if (lJe > 1) { 
    u_matrix(n,1,r[1],r[2],mat_c);
    for (j=1; j<=lJe-1; j++) {
      
      u_matrix(n,j,r[j],r[j+1],mat_);

      copy_2x2(mat_d,mat_c);
      matrix_prod_2x2(mat_d,mat_,mat_c);
    }
    u_matrix(n,lJe,r[lJe],lr_e,mat_a);

    matrix_prod_2x2(mat_c,mat_a,mat_d);
    Rpr=deep?
      mat_d[2][1]:mat_d[2][2];
  } else { 
    if (!lJe){

      Rpr=1.0f; 
    }
    else {
      u_matrix(n,1,r[1],lr_e,mat_a);
      Rpr= deep?
	mat_a[2][1]:mat_a[2][2];/* [2][1] */
    }    
  }
  if (n > rbsize) {
    printf("-"); fflush(NULL);
    rbsize +=100;
    Rbottom = (double*)realloc(Rbottom,sizeof(double)*(rbsize+1));
    Rright = (double*)realloc(Rright,sizeof(double)*(rbsize+1));
    legendre = (double*)realloc(legendre,sizeof(double)*(rbsize+2));
    for (j=rbsize-100; j <= rbsize; j++)
      Rbottom[j]=0.0f; 

  }

  /* second of r_e and r_o through to the center. */
  /* can be automated only for the shallow case. */
  if (deep || (Rright[n]==0.0f)) {  
    u_matrix(n,lJo,lr_o,r[lJo+1],mat_c);

    for (j=lJo+1; j<=layers; j++) {
      copy_2x2(mat_d,mat_c);
      u_matrix(n,j,r[j],r[j+1],mat_);
      matrix_prod_2x2(mat_d,mat_,mat_c);
    }
    if (!deep)
      Rright[n] =mat_c[2][2];

    Rpr*=  deep?mat_c[1][2]:mat_c[2][2];/* for deep 1 2 */
  } else 
    Rpr*=Rright[n];

  if (Rbottom[n] == 0.0f) {  
    u_matrix(n,1,r[1],r[2],mat_c);
    for (j=2; j<=layers; j++) {
      copy_2x2(mat_d,mat_c);
      u_matrix(n,j,r[j],r[j+1],mat_);
      matrix_prod_2x2(mat_d,mat_,mat_c);
    }
    Rbottom[n] = -mat_c[2][2];

  }
  Rpr/=lr_e*lr_e*Rbottom[n]*eps[Jo];


  Rpr*=pow(lr_o/r[lJo],nu_lJo)/pow(lr_e/r[lJe],nu_lJe);       
  for (j=lJe; j<=lJo-1; j++) {
    nu_=nu(j,n);
    Rpr *=lambda(j,nu_);
  }

  return Rpr;

}


int find_Jo(double r_o) { /* see 27*/ 
  int temp=1;
  while (r[temp]>=r_o)
    temp++;
  temp--;

  return temp;
}

int find_Je(double r_e) {/* see 27*/ 
  int temp=1;
  while (r[temp]>r_e)
    temp++;
  temp--;

  return temp;
}

#if 1

double Smono(double costheta, 
	     double r_o, double r_e) { 
  /* equation 4 - monopole! */
  double S,Rn,Sprev; 
  int n,flag;
  
  S=Sprev=0.0f;
  if (r_o<0.0f) 
    return 0.0f;
  legendre[0]=1;
  legendre[1]=costheta; 
  for (n=1,flag=1;flag;n++){
    Sprev=S;
    Rn=R(n,r_o,r_e,find_Jo(r_o),find_Je(r_e));
    if (n>1)
      legendre[n] =
	((2*n-1)*costheta*legendre[n-1]-
	 (n-1)*legendre[n-2])/n;    
    S+=(2*n+1)*Rn*legendre[n];
    if (n>minimum_terms &&(fabs((S-Sprev)/S) < tolerance))
      flag = 0;
    if (!finite(S)) {

      return Sprev; 
    }
  }
  return S;

  

}

double S0 (double costheta,
	   double r_o, double r_e) { 
  /* equation 57   */
  /* tangential dipoles. */
  double S,Ra,P,F0,Lambda,A,Rn;
  double lp, Sprev;
  int i,n,flag;
  int Jo,Je;
  Jo=find_Jo(r_o);  Je=find_Je(r_e);
  if (!slow && (r_e>r_o) ){

    if (isotropic) { 
      Lambda = r_o/r_e; 
    }else { 
      Lambda =
	pow(r_o/r[Jo],alpha(Jo))/
	pow(r_e/r[Je],alpha(Je));
      A =
	pow(r_o/r[Jo],0.5f*(alpha(Jo)-1))/
	pow(r_e/r[Je],0.5f*(alpha(Je)-1)); 
      for (i=Jo-1; i>= Je; i-- ) {
	
	lp = lambda(i,1.0f);
	Lambda *= pow(lp, alpha(i));
	A *= pow(lp,
		 0.5f*(alpha(i)-1)
		 );
      }
    }

    F0 = -A/(r_e*beta(Jo));
    if (Je != Jo)
      for (i=Je ; i<= Jo-1; i++ ) 
	F0 *= 2*beta(i+1)/(beta(i) + beta(i+1));
    
    
    Ra = 1-2*Lambda*costheta + Lambda*Lambda;
    lp = 1.0f;
  }
  S = (slow ||(r_e<r_o))?0.0f: F0*Lambda/(r_o*Ra*sqrt(Ra));
  legendre[0]=1;
  legendre[1]=costheta; 

  for (n = flag = 1; flag > 0 ; n++) {
    Sprev = S;
    if (n>1)
      legendre[n] =
	((2*n-1)*costheta*legendre[n-1]-
	 (n-1)*legendre[n-2])/n;
    /* this is how it gets differentiated */
    P = (n*costheta*legendre[n]-n*legendre[n-1])/(costheta*costheta-1);
    Rn = R(n,r_o,r_e,find_Jo(r_o),find_Je(r_e));

    if (slow ||(r_e<r_o)) 
      S+= (2*n+1)*Rn*P/r_o;
    else{
      lp *= Lambda;
      S += ((2*n + 1) * Rn - F0*lp)*P/r_o;
    }
    if (n>minimum_terms &&(fabs((S-Sprev)/S) < tolerance))
      flag = 0;
    if (isnan(S)) {
     
      return Sprev; 
    }
  }
  return S;
}


double S1 (double costheta,double r_o, double r_e) { 
  /* equation 58 */  
  /* radial dipoles */
  double S,F0,F1,Ra,P,Lambda,A,Rn;
  double lp, Sprev;

  int i,n,flag;
  int Jo,Je;
  Jo=find_Jo(r_o);  Je=find_Je(r_e);

  if (!slow && (r_e>r_o)){
    if (isotropic) { 
      Lambda = r_o/r_e; 
      A=1.0f;
    }else { 
      Lambda =
	pow(r_o/r[Jo],alpha(Jo))/
	pow(r_e/r[Je],alpha(Je));
      A =
	pow(r_o/r[Jo],0.5f*(alpha(Jo)-1))/
	pow(r_e/r[Je],0.5f*(alpha(Je)-1)); 
      for (i=Jo-1; i>= Je; i-- ) {
	
	Lambda *= lambda(i, alpha(i));
	A *= lambda(i,
		    0.5f*(alpha(i)-1)
		    );
      }
    }
    F0 = -A/(r_e*beta(Jo));
    for (i=Je ; i<= Jo-1; i++ ) {
      F0 *= 2*beta(i+1)/(beta(i) + beta(i+1));
    }
    F1 =
      F0 *
      sqrt(eps[Jo]*eta[Jo]) /(r_o*eps[Jo]);
    
    Ra = 1.0f -2*Lambda*costheta + Lambda*Lambda;
    Ra *=sqrt(Ra);  
    lp = 1.0f;
  }
  S = (slow ||(r_e<r_o)) ?0.0f:F1*Lambda*(costheta-Lambda) /Ra;

  legendre[0]=1;
  legendre[1]=costheta; 
  for (n = 1; ; n++) {
    Sprev = S;

    Rn = Rprime(n,r_o,r_e,Jo,Je);

    if (n>1)
      legendre[n] =
	((2*n-1)*costheta*legendre[n-1]-
	 (n-1)*legendre[n-2])/n;
    if (slow ||(r_e<r_o)) 
      S+= (2*n+1)*Rn*legendre[n];
    else { 
      lp *= Lambda;
      S += ((2*n + 1) * Rn - F1*n*lp)*legendre[n];
    }
    if (isnan(S))
      return Sprev; 
    if (n>minimum_terms && (fabs((S-Sprev)/S) < tolerance))
      break;     
  }
  return S;
}
double psi_function(double r_o, double r_e,
		 double costheta, double cosphi,
		 double Mrad1, double Mtan1){

  double psi=0.0f;
  double Sa,Sb; 
  /* this worked for tangentials on 
     June 1st, with S0 doing tangentials. */
  find_Je(r_e);
  if (mono>0.0f) {

    if (Mtan1 != 0.0f){ 
      fprintf(stderr,
	      "Can't do tangential monopole pairs yet.\n");
      exit(-22);
    }
    if (Mrad1 != 0.0f){
      Sa=Smono(costheta,r_o,r_e);
#if 0
      Sb=Smono(costheta,r_o-mono*dz,r_e);
#else 
      Sb=0.0f;
#endif
      psi+= Mrad1*(Sa-Sb);

    }
  } else {
    if (Mtan1 != 0.0f) 
      psi += Mtan1* cosphi* S0(costheta,r_o,r_e);
    if (Mrad1 != 0.0f)
      psi += Mrad1* S1 (costheta,r_o,r_e);
  }
  return psi/FOUR_PI;
}
#endif
double r_e,r_o;
#define MAXLINE 1024
void calculate_sensors(double mrad ,double mtan,
		       char *infilename){
  FILE * infile;
  char line[MAXLINE];
  int numsensors=0,n=0; 
  double x,y,z,az,el,vi;
  infile = fopen(infilename,"r"); 
  if (!infile) {
    fprintf(stderr,"%s not found.\n",infilename); 
    exit(-1); 
  }
  rbsize = 1000; 
  r_o=r[0];
  find_Jo(r_o);
  Rbottom = (double *)malloc(sizeof(double)*(rbsize+1)); 
  Rright = (double *)malloc(sizeof(double)*(rbsize+1)); 
  legendre = (double *)malloc(sizeof(double)*(rbsize+2)); 
  for (n=0; n <= rbsize; n++)
    Rbottom[n]=0.0f; 

  if (filemode>1) { 
    if (!fscanf(infile,"%d\n",&numsensors))
      fprintf(stderr,"File lacks correct format.\n");
    for (n=0;n<numsensors;n++) {
      if (!fscanf(infile,"%lf %lf %lf %*f %*f %*f\n",&x,&y,&z))
	fprintf(stderr,"File goof.\n");; 
      /* using the dx,dy,dz variables helps 
	 with unit conversions. */
      r_e = sqrt(x*x*dx*dx+y*y*dy*dy+z*z*dz*dz); 
      el=dz*z/r_e;    
      az=dy*y/r_e;
      vi = psi_function(r_o,r_e,el,az,mrad,mtan); 
      printf("%f\n",vi);
    }
  } else { 
    /* if it's a PPI file: */

    int ncom,c,i,j,ncol,laten,N_sensors;
    float ftemp;
    /* skip first line*/
    fgets(line,MAXLINE,infile);
    
    /* read number of lines in the header */
    fgets(line,MAXLINE,infile);
    sscanf(line,"%4d%c",&ncom,&c);
    
    /* look for latencies and noise in header */
    laten = 0; /* laten is set to 1 if the next input line 
		  will contain the latency values */
    for (j=0; j<ncom ; j++)
      {
	fgets(line,MAXLINE,infile);
      }
    
    /* read number of columns of data */
    fgets(line,MAXLINE,infile);
    sscanf(line,"%4d",&ncol);
    /* measure page size */
    N_sensors = 0;

    while (fgets(line,MAXLINE,infile))
      if (strlen(line) > (size_t)10) (N_sensors)++;

    rewind(infile);

    /* Now skip the header, and read the sensor and amplitude data. */
    for (i=0; i<ncom+3; i++)
      fgets(line,MAXLINE,infile);
    
    for (i=0; i<N_sensors; i++) {
      fscanf(infile, "%*f%*f%*f%*f%*f%*f%");

      fscanf(infile, "%lf %lf %lf", &x,&y,&z);
      x+=sx;
      y+=sy;
      z+=sz;
      for (j=0; j<3; j++)
	{
	  fscanf(infile, "%*f");
	}
      for (j=0; j<ncol; j++)
	{
	  fscanf(infile, "%*f",&ftemp);
	}

      r_e = sqrt(x*x*dx*dx+y*y*dy*dy+z*z*dz*dz); 
      el=z*dz/r_e;    
      az=y*dy/r_e;
      vi = psi_function(r_o,r_e,el,az,mrad,mtan); 
      printf("%f\n",vi);
    }
    
  } 
  
  fclose(infile);
}


void record_volume(
                    double mrad ,double mtan,
                    char *outfilename)
{
  double xc,yc,zc,az,el;
  int kb,dummy;
  double r_oe;
  float temp_float; 
  FILE *surfile;
  surfile= fopen( outfilename,"w");
  printf("Recording the surface.\n");
  rbsize = 1000; 
  Rbottom = (double *)malloc(sizeof(double)*(rbsize+1)); 
  Rright = (double *)malloc(sizeof(double)*(rbsize+1)); 
  legendre = (double *)malloc(sizeof(double)*(rbsize+2)); 
  for (dummy=0; dummy <= rbsize; dummy++)
    Rbottom[dummy]=0.0f; 

  r_o=r[0];
  r[0]=-10000.0f;/* trap! */
  find_Jo(r_o);
  dummy=0;
  for (curr_i=start_i,curr_j=start_j,curr_k=start_k;
       curr_k<=K;
       curr_k++) {
    zc=dz*((double)curr_k-(double)(K/2)-sz);
    printf("%d",curr_k%10);
    if (dummy>0)
      printf(".");
    fflush(NULL);

    for (;curr_j<=J;curr_j++){

      yc=dy*((double)curr_j-(double)(J/2)-sy);
      for (;curr_i<=I;curr_i++){
	xc=dx*((double)curr_i-(double)(I/2)-sx);
	
	r_e=sqrt(xc*xc+yc*yc+zc*zc);	
	//	r_oe = sqrt(xc*xc+yc*yc+(zc-r_o)*(zc-r_o));
	/* we're generating our own volumes now */
	if (/*(voxel_radius >r_o)&&*/
	    (r_e<= r[1])) {
	  dummy++;
	  /* no orientation switches.
	     We be standardising for the bakeoff.
	     Dipole is always on the slowest changing axis,
	     and pointing to the second slowest in the case 
	     of tangentials. */

	  el=zc/r_e;
	      
	  az=yc/r_e;

	  temp_float=(float)psi_function(r_o,r_e,el,az,mrad,mtan);

	} else 
	  temp_float = 0.0f; 
	if (isnan(temp_float))
	  temp_float = 0.0f; 
	fwrite(&temp_float, sizeof(float) , 1, surfile);

      }
      curr_i=1;
    }
    curr_j=1;
  }
  
  fclose(surfile);
  printf("\n");
   
}
extern char *optarg;
extern int optind, opterr, optopt;

int main(int argc, char *argv[])
{
  int i,j,k;
  int surfmode,index,option_index;
  FILE *infile;

  char * comma;
  char c;

  struct option long_options[] =
  {
    {"I",1,0,'I'},
    {"J",1,0,'J'},
    {"K",1,0,'K'},
    {"filemode",1,0,'f'},
    {"sensorlist",0,0,'F'},
    {"or",1,0,'o'},
    {"orientation",1,0,'o'},
    {"resolution",1,0,'d'},
    {"recenter",1,0,'S'},
    {"dxdydz",1,0,'d'},
    {"layers",1,0,'l'},
    {"eta",1,0,'h'},
    {"eps",1,0,'e'},
    {"radii",1,0,'r'},
    {"Mrad",1,0,'R'},  
    {"Mtan",1,0,'T'},
    {"minimum",1,0,'m'},
    {"Tolerance",1,0,'t'},
    {"i",1,0,'i'},
    {"j",1,0,'j'},
    {"k",1,0,'k'},
    {0, 0, 0, 0}
  };

  signal(SIGHUP,oopsie); 
  signal(SIGINT,oopsie); 

  layers=4;

  I=100;
  
  J=100;
  K=100;

  r=fvector(layers+1);
  eps=fvector(layers);
  eta=fvector(layers);



  Mrad=1000;
  Mtan=0;
  filemode=1;

  /* Because I am an MKS bigot, the units are as follow: */
  while(1) {
    c = getopt_long (argc, argv, "I:J:K:R:T:r:e:h:f:o:l:d:t:m:S:scFH",
		     long_options, &option_index);
    if (c==-1) 
      break; 
    switch(c){
    case 'H':
      fprintf(stderr,"Usage:\n"
	      "demunck filename -I I -J J -K K \n"
	      "-l (number of layers) -r [radii]\n"
	      "-e [radial conductivity] -h [tangential]\n"
	      "-R (radial dipole) -T (tangential) \n"
	      "-d [resolution] -S [displacement] \n"
	      "Bunch of others.\n");
      exit(0);
    case 'm':
      minimum_terms = atoi(optarg);
      break;
    case 'c': 
      singleslice=1;
      break;
    case 's': /* forget the addition-subtraction speedup. */
      slow =1;
      break;
    case 'd': /* resolution of our cube in meters */
      sscanf(optarg,"%lf,%lf,%lf",&dx,&dy,&dz);
      break;
    case 'f': 
      filemode=atoi(optarg);
      break;
    case 'I':   /* number of voxels (for when we generate cubes) */
      I=atoi(optarg);
      break; 
    case 'J':
      J=atoi(optarg);
      break;
    case 'K':
      K=atoi(optarg);
      break; 
    case 'i': /* starting dimensions */
      start_i=atoi(optarg);
      break; 
    case 'j':
      start_j=atoi(optarg);
      break;
    case 'k':
      start_k=atoi(optarg);
      break; 
    case 'o':
      mono=atof(optarg);
      break;
    case 'S':
      sscanf(optarg,"%lf,%lf,%lf",&sx,&sy,&sz);
      break;

    case 'l': /* number of spherical layers */
        layers=atoi(optarg);
	r=re_fvector(r,layers+1);
	eps=re_fvector(eps,layers);
	eta=re_fvector(eta,layers);
	break;
    case 'F':
      sensorlist =1; 
      break;
    case 'R':
      Mrad=atof(optarg);/* ampere*meters */
      break;
    case 't': 
      tolerance = atof(optarg);
      break; 
    case 'T':
      Mtan=atof(optarg); /* ampere*meters */
      break;
    case 'r': /* radii, r0 being radius of dipole, r1-N 
		 radii of the layers, descending order */
      index=0; 
      comma=optarg;
      while ((index <=layers) && (comma != NULL)) {
	sscanf(comma,"%lf",&r[index]); /* meters */
	comma = (char*)strchr(comma,',');
	if (comma!=NULL) 
	  comma++; 
	index++;
      }
      r[layers+1]=0.0f;
      break;
    case 'h': /* tangential conductivity (S/m) */
      index=1; 
      comma=optarg;
      while (comma != NULL) {
	sscanf(comma,"%lf",&eta[index]); /* siemens per meter */
	comma = (char*)strchr(comma,',');
	if (comma!=NULL) 
	  comma++; 

	index++;
      }

      break;
    case 'e':/* radial conductivity (S/m) */
      index=1; 
      comma=optarg;
      while (comma != NULL) {
	sscanf(comma,"%lf",&eps[index]);/* siemens per meter */
	comma =(char*) strchr(comma,',');
	if (comma!=NULL) 
	  comma++; 
	index++;
      }

    }
  }
  fprintf(stderr,"I %d. ",I);

  fprintf(stderr,"J %d. ",J);

  fprintf(stderr,"K %d. ",K);

  fprintf(stderr," filemode %d\n",filemode);
  fflush(NULL);
  fprintf(stderr,"Mrad: %f Mtan: %f \n",
         Mrad, Mtan);
  for (index=1; index<=layers; index++) { 
    if (eps[index] != eta[index])
      isotropic =0 ;
  }
  fprintf(stderr,"Resolution: %f %f %f\n",dx,dy,dz);
  fprintf(stderr,"radii: ");
  for (i=0; i<=layers+1; i++)
    fprintf(stderr,"%f ",
	   r[i]);
  fprintf(stderr,"\neta: ");
  for (i=1; i<=layers; i++)
    fprintf(stderr,"%f ",
	   eta[i]);         
  fprintf(stderr,"\neps: ");
  for (i=1; i<=layers; i++)
    fprintf(stderr,"%f ",
	   eps[i]);         
  fprintf(stderr,"\n");
  fprintf(stderr,"S %f %f %f\n", sx,sy,sz);
  if (tolerance == 0.0f) {
    fprintf(stderr, "Nonzero tolerance value req'd. \n");
    exit(-1);
  }

  if (!sensorlist)
    record_volume(Mrad,Mtan,argv[argc-1]);
  else 
    calculate_sensors(Mrad,Mtan,argv[argc-1]);

  return 0; 
}
/*fid=fopen('/var/tmp/homo2.flt','r','ieee-le'); a=fread(fid,inf,'float'); fclose(fid); a=reshape(a,[150 150 150]); 
 pcolor(atan(squeeze(a(:,75,:)))); shading flat

*/