#define EXTERN extern  /*when read def.h make all EXTERN variables extern*/
#include "global.h"  // global variables  color functions  (external variable)

#define TINY 1.0e-20


float sqr(float x);
void orthtofr(struct griddata& den1,float xx[3],float fr[3]); 
void frtoorth(struct griddata& den1,float fr[3],float orth[3]); 
void xtof(struct griddata& den1,float xx[3],float ijk[3]); /*convert from world co-ord to floating i,j,k co-ord*/
void xtoi(struct griddata& den1,float xx[3],int ijk[3]); /*convert from world co-ord to interger i,j,k co-ord*/
void imcv(struct griddata& den1,int xi[3],float sx[3]);  /*convert from xi,xj (in fractional floating point)
	 to world co-ord, identical to cv but uses vectors*/
void mcv(struct griddata& den1,float xi[3],float sx[3]);  /*convert from xi,xj (in fractional floating point)
	 to world co-ord, identical to cv but uses vectors*/
void imat2(int** &mat,int n1,int n2);// allocates array 2 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
void freeimat2(int** &mat,int n1,int n2);// free memory array 3 dim array (mat)
void freefmat2(float** &mat,int n1,int n2);// free memory array 3 dim array (mat)
void fmat2(float** &mat,int n1,int n2);// allocates array 2 dim array (mat)
         // of type float - n1,n2,n3,n4 size of each dimension
void cmat2(char** &mat,int n1,int n2);// allocates array 2 dim array (mat)
         // of type char - n1,n2,n3,n4 size of each dimension
void cmat3(char*** &mat,int n1,int n2,int n3);// allocates array 3 dim array (mat)
         // of type char - n1,n2,n3,n4 size of each dimension
void imat3(int*** &mat,int n1,int n2,int n3);// allocates array 3 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
void freeimat3(int*** &mat,int n1,int n2,int n3);// free memory array 3 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
        // should free in reverse order of allocating
void fmat3(float*** &mat,int n1,int n2,int n3);// allocates array 3 dim array (mat)
         // of type  signed char - n1,n2,n3,n4 size of each dimension
void freefmat3(float*** &mat,int n1,int n2,int n3);// free memory array 3 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
        // should free in reverse order of allocating
void simat3(short int*** &mat,int n1,int n2,int n3);// allocates array 3 dim array (mat)
         // of type short int - n1,n2,n3,n4 size of each dimension
void scmat3(signed char*** &mat,int n1,int n2,int n3);// allocates array 3 dim array (mat)
         // of type  signed char - n1,n2,n3,n4 size of each dimension
void unscmat3(unsigned char*** &mat,int n1,int n2,int n3);// allocates array 3 dim array (mat)
         // of type unsigned signed char - n1,n2,n3,n4 size of each dimension
void freeunscmat3(unsigned char*** &mat,int n1,int n2,int n3);// free memory array 3 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
        // should free in reverse order of allocating
void freescmat3(signed char*** &mat,int n1,int n2,int n3);// free memory array 3 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
        // should free in reverse order of allocating
void scmat4(signed char**** &mat,int n1,int n2,int n3,int n4);// allocates array 4 dim array (mat)
         // of type signed char - n1,n2,n3,n4 size of each dimension
void freescmat4(signed char**** &mat,int n1,int n2,int n3);// free memory array 4 dim array (mat)
void fmat4(float**** &mat,int n1,int n2,int n3,int n4);// allocates array 4 dim array (mat)
         // of type float - n1,n2,n3,n4 size of each dimension
void freefmat4(float**** &mat,int n1,int n2,int n3,int n4);// free memory array 4 dim array (mat)
void fmat5(float***** &mat,int n1,int n2,int n3,int n4,int n5);// allocates array 5 dim array (mat)
         // of type float - n1,n2,n3,n4 n5 size of each dimension
void simat5(short int***** &mat,int n1,int n2,int n3,int n4,int n5);// allocates array 5 dim array (mat)
         // of type short - n1,n2,n3,n4 n5 size of each dimension
int d3tod1(int cmax[],int i,int j,int k); // convert a 3 dimsional positiion to a 1 dimensional one
void d1tod3(int cmax[],int d1,int d3v[3]);
void d1tow(struct griddata& den1,int cmax[],int ipos,float fvec[3]);
void planevec(float cn[3],float nca[3],float cdn,float sn,
                float r,float th,float th1,float vorigin[3],float vec[3]);
float dot(float vec1[],float vec2[]); /*returns v1 dot v2*/
void crossprod(float va[3],float vb[3],float vcross[3]); /* (float) cross prod. 
         of 2 float  vectors*/
void rotaxis(int de[3][3],int ep[3][3][3],float phi,float nca[3],
                       float rotm[3][3]);
void tranm(float v[3],float t[3]); /*return vector v=v-t*/
void trannewm(float v[3],float vnew[3],float t[3]); /*return vector vnew=v-t*/
void tran(float v[3],float t[3]); /*return vector v=v+t*/
void trannew(float v[3],float vnew[3],float t[3]); /*return vector vnew=v+t*/
void multvec(float m[3][3],float v1[],float v[]); /* returns vector v = m*v1*/
void imultvec(int m[3][3],float v1[],float v[]); /* returns vector v = m*v1*/
float distv(float vec1[],float vec2[]); // Returns the distance between vec 1 and vec 2
float determ(float a[3][3]);
void nrerror(const char error_text[]);
void findchivec(float u[3],float v[3],float w1[3],float *chi);
void findchi(float p0[3],float p1[3],float p2[3],float p3[3],float *chi);
void rotateset(float orig[3],float rpt[3],float rom[3][3]);
void rotatesetnew(float orig[3],float rpt[3],float tempv[3],float rom[3][3]);
void rotatept(struct geometry geo,float rvec[3],float orig[3],float ang,
      float rpt[3]);
void rottranvec(float rot[3][3],float tran[3],float v1[3],float v2[3]);
void printmat3(char aname[20],float a[3][3]);
void printmat6(char aname[20],float a[6][6]);
void mult11_61(float aa,float bb[6],float cc[6]);
void tran66(float aa[6][6],float cc[6][6]);
void add66_66(float aa[6][6],float bb[6][6],float cc[6][6]);
void mult33_31(float aa[3][3],float bb[3],float cc[3]);
void mult33_33(float aa[3][3],float bb[3][3],float cc[3][3]);
void mult66_66(float aa[6][6],float bb[6][6],float cc[6][6]);
void mult66_61(float aa[6][6],float bb[6],float cc[6]);
void mult16_66(float aa[6],float bb[6][6],float cc[6]);
void mult16_61(float aa[6],float bb[6],float *cc);
void mult61_16(float aa[6],float bb[6],float cc[6][6]);
void mymatinv(float a[6][6], float ai[6][6]);
void mymatinv3(float a[3][3], float ai[3][3]);
void myludcmp(float a[6][6],int indx[6],int *d);
void mylubksb(float a[6][6],int indx[6],float b[]);
void solvelineq(int ndim,float a[5][5],float b[]);
void itoa(int ii,char ch[]);
void rankqual(int maxrank,int ranki[],int ranki2[],float rankf[],int ivalue,float qual,
     int gapnum);
void rankqual3(int maxrank,int ranki[],int ranki2[],int ranki3[],float rankf[],int ivalue,
     float qual,int gapnum1,int gapnum2);
void newrankqual(float qual,int ranki[4][MAXRANK],int intprop[4],
        float rankf[MAXRANK]);
void orderlist(int del[],int order[]);

 

 

