/*                ***********NOTICE**********                              */
/*                                                                         */
/*  The EMAP finite element modeling codes were created at the             */
/*  University of Missouri-Rolla Electromagnetic Compatibility Laboratory. */
/*  They are intended for use in teaching or research.  They may be freely */
/*  copied and distributed PROVIDED THE CODE IS NOT MODIFIED IN ANY MANNER.*/
/*                                                                         */
/*  Suggested modifications or questions about these codes can be          */
/*  directed to Dr. Todd Hubing, Department of Electrical Engineering      */
/*  University of Missouri-Rolla, Rolla, MO  65401.  Principal authors     */
/*  of these codes are Mr. Mohammad Ali and Mr. Girish Bhat of the         */
/*  University of Missouri-Rolla.                                          */
/*                                                                         */
/*  Neither the authors nor the University of Missouri makes any warranty, */
/*  express or implied, or assumes any legal responsibility for the        */
/*  accuracy, completeness or usefulness of these codes or any information */
/*  distributed with these codes.                                          */
/*                                                                         */
/*  2/2/95                                                                 */

/***************************************************************************
*   PROGRAM     :  Emap2.c (Version 2.02)
*   LAST UPDATE :  February 15, 1995
*   DESCRIPTION :  A 3-D Finite Element Modeling Code for Analyzing
                   Time Varying Complex Electromagnetic fields.
                   (a scalar code with spurious mode suppression)
*   INPUT FILE  :  emap2.in or or 1st argument of emap command
*   OUTPUT FILE :  file(s) specified by the input file
*
****************************************************************************/

/**************************** Include Files ********************************/

#include<stdio.h>
#include<signal.h>
#include<string.h>
#include<stdlib.h>
#include<sys/time.h>
#include<math.h>

/******************** Dynamic Memory Allocation Functions ******************/

int *INT_Vector(nh)
/* Allocates an int vector with range [0 ... nh-1]  */
int nh;
{
int *v;
v=(int *)malloc((unsigned) nh*sizeof(int));
return v;
}
double *FLOAT_Vector(nh)
/* Allocates a double vector with range [0 ... nh-1] */
int nh;
{
double *v;
v=(double *)malloc((unsigned) nh*sizeof(double));
return v;
}

double **FLOAT_Matrix(row,column)
/*Allocates a double matrix with range [0..rw-1][0..cl-1]*/
int row,column;
{
int i;double **m;
m=(double **) malloc((unsigned) row*sizeof(double*));
for(i=0;i<row;i++)
{m[i]=(double *) malloc((unsigned) column*sizeof(double));}
return m;
}
int **INT_Matrix(row,column)
/*Allocates an int matrix with range [0..rw-1][0..cl-1]*/
int row,column;
{
int i;int **m;
m=(int **) malloc((unsigned) row*sizeof(int));
for(i=0;i<row;i++)
{m[i]=(int *) malloc((unsigned) column*sizeof(int));}
return m;
}

/************************* Default Parameters ******************************/

#define  GBL_Matrix_MaxColNos 43
#define  FRCD_Matrix_MaxColNos 10

/*********************** Function Prototypes *******************************/

void ParameterInfo();
void Assign_Global_Coordinates();
void Assign_hexahedron_NodeNumbering();
void Find_CoFactor();
void Find_Volume_Determinant();
void HexaHedra_SubDivision_1();
void HexaHedra_SubDivision_2();
void TetraHedron_SubMatrix();
void HexaHedron_Matrix();
void HexaHedron_link_1();
void HexaHedron_link_2();
void Find_Global_Matrix();
void Read_Input_Pass_1();
void Read_Input_Pass_2();
void Fill_Vectors();
void Partition_Global_Matrix();
void Find_RHS_Vector();
void Band_Matrix_Solver();
void Produce_Output();
void EfieldatNode();
void Count_Half_BandWidth();
void swap();
/***************************************************************************/

double Determinant_TetraMatrix();
int Locate_NZ_Element_Column();
int Srch_NZ_Element_Column();

/*********************** Fixed Matrix Pointers *****************************/

double cofs[4][4];
double Volume_Cord[4][4],VVolume_Cord[4][4];
int   tetrnod[5][4],tetlink[5][4],connect1[8][6];
double Tet_Matrix_Element[5][12][12],Hex_Matrix_Element[24][24];

/************************** Global Pointers ********************************/

double *REL_Permittivity,*FORC_Val,**Dat_Buffer,**GBL_Matrix_Data;
double **MeshNodeCord,**FORC_Matrix_Data,*RHS_Vector,**HalfBand_Matrix;
int **GBL_Matrix_ColNos,*GBL_Matrix_Index;
int **FORC_Matrix_ColNos,*UN_Vector,*FORC_Vector,**HexHed;
int *Loc_UN_Node,*Loc_UN_NodeCord;
FILE *InF, *OutF_0, *OutF_1, *OutF_2, *OutF_3;

/************************** Global Variables *******************************/

int Dimension_X,Dimension_Y,Dimension_Z;
double OperateFreq,CellDimension;
double DivisorX,DivisorY,DivisorZ;
int Min_X=9999,Min_Y=9999,Min_Z=9999,Max_X=-9999,Max_Y=-9999,Max_Z=-9999;
double det=0.0, WaveNumber, freq;
int t, cn=0;
int TotNum_GBL_Nodes,TotNum_HexaHedron;
struct timeval;
int Half_BandWidth=0,TotNum_Unknown_Var=0,TotNum_Known_Var=0;
char Out_FileName0[20],Out_FileName1[20],Out_FileName2[20],Out_FileName3[20];

/************************** Main Function **********************************/

main(argc, argv)

int argc;
char *argv[];

{
int    i,j,GBL_Matrix_RowNos;
char   Ifile[20];
        if (argc > 2) {
        fprintf(stderr,"Usage: Emap2 [input_file] \n"); exit(1);}
        if (argc < 2) {  argv[1] = "emap2.in";   strcpy(Ifile, argv[1]); }
        else                                    strcpy(Ifile, argv[1]);
        fprintf(stdout, "\nEMAP2 Version 2.02\n\n");
        fprintf(stdout, "The input will be read from the file:  %s\n", Ifile);
        InF = fopen(argv[1], "r");
        if (InF == NULL) {
        fprintf(stderr, " Error: Can't open input file.\n"); exit(1); }
        
Read_Input_Pass_1();
ParameterInfo();

/* Allocates Dynamic memory */
HexHed=INT_Matrix(TotNum_HexaHedron,8);
MeshNodeCord=FLOAT_Matrix(TotNum_GBL_Nodes,3);
REL_Permittivity=FLOAT_Vector(TotNum_GBL_Nodes);
Dat_Buffer=FLOAT_Matrix(TotNum_GBL_Nodes,3);

Read_Input_Pass_2();
fprintf(stderr,"TotNum_Unknown_Var=%d\n", TotNum_Unknown_Var);
fprintf(stderr,"TotNum_Known_Var=%d\n", TotNum_Known_Var);

Assign_Global_Coordinates();
Assign_hexahedron_NodeNumbering();

/* Allocates Dynamic memory */
 UN_Vector=INT_Vector(TotNum_Unknown_Var);
 FORC_Vector=INT_Vector(TotNum_Known_Var);
 FORC_Val=FLOAT_Vector(TotNum_Known_Var);
 Loc_UN_Node=INT_Vector(TotNum_Unknown_Var);
 Loc_UN_NodeCord=INT_Vector(TotNum_Unknown_Var);
 FORC_Matrix_Data=FLOAT_Matrix(TotNum_Unknown_Var,FRCD_Matrix_MaxColNos);
 FORC_Matrix_ColNos=INT_Matrix(TotNum_Unknown_Var,FRCD_Matrix_MaxColNos);
 RHS_Vector=FLOAT_Vector(TotNum_Unknown_Var);
 GBL_Matrix_RowNos=TotNum_GBL_Nodes*3;
 GBL_Matrix_Data=FLOAT_Matrix(GBL_Matrix_RowNos,GBL_Matrix_MaxColNos);
 GBL_Matrix_ColNos=INT_Matrix(GBL_Matrix_RowNos,GBL_Matrix_MaxColNos);
 GBL_Matrix_Index=INT_Vector(GBL_Matrix_RowNos);


/*  Initialize Dynamically Allocated Variables  */

for(i=0;i<=3*TotNum_GBL_Nodes-1;i++) {
  GBL_Matrix_Index[i]=0;
  for(j=0;j<=GBL_Matrix_MaxColNos-1;j++) {
    GBL_Matrix_Data[i][j]=0.0;
    GBL_Matrix_ColNos[i][j]=0; }}
   for(i=0;i<=23;i++) for(j=0;j<=23;j++)
   { Hex_Matrix_Element[i][j]=0.0;}
   for(i=0;i<=11;i++)
    for(j=0;j<=11;j++)
      for(t=0;t<=4;t++)
       { Tet_Matrix_Element[t][i][j]=0.0;}

Fill_Vectors();

for(cn=0;cn<=(TotNum_HexaHedron-1);cn++) {      
 if(cn%2 ==1) HexaHedra_SubDivision_2();
 else         HexaHedra_SubDivision_1();

   for(i=0;i<=11;i++)
    for(j=0;j<=11;j++)
      for(t=0;t<=4;t++)
       { Tet_Matrix_Element[t][i][j]=0.0;}
  for(t=0;t<=4;t++)          {
    Find_Volume_Determinant();         
    Find_CoFactor();          
    TetraHedron_SubMatrix(); }

    HexaHedron_Matrix(); 

   if(cn%2 ==1) HexaHedron_link_2();
   else         HexaHedron_link_1();
   Find_Global_Matrix();

    } /* end of cn loop */

Count_Half_BandWidth();

/* Allocates Dynamic memory */
HalfBand_Matrix=FLOAT_Matrix(TotNum_Unknown_Var+1,Half_BandWidth+1);

for(i=1;i<=TotNum_Unknown_Var;i++)
 for(j=1;j<=Half_BandWidth;j++)
 {HalfBand_Matrix[i][j] =0.0;}

Partition_Global_Matrix();
Find_RHS_Vector();

Band_Matrix_Solver();
Produce_Output();
return(0);
}

/*********************End of Main Program **********************************/
/***************************************************************************
*
*   FUNCTION        : ParameterInfo()
*   DESCRIPTION     : Initializes global variables. 
*
****************************************************************************/
void ParameterInfo()
{
double FreeSpaceVel,AbsPermeable,AbsPermitt,WaveLength, UCellLen;
 
 AbsPermeable  = 1.25663706144E-06;
 AbsPermitt    = 8.8542E-12;
 FreeSpaceVel  = 1.0/sqrt(AbsPermeable*AbsPermitt);
 WaveLength = FreeSpaceVel/(OperateFreq*1.0E+06);
 WaveNumber = 2.0*M_PI/WaveLength;

 Dimension_X=Max_X-Min_X; Dimension_Y=Max_Y-Min_Y; Dimension_Z=Max_Z-Min_Z;
 TotNum_HexaHedron=Dimension_X*Dimension_Y*Dimension_Z;
 TotNum_GBL_Nodes=(Dimension_X+1)*(Dimension_Y+1)*(Dimension_Z+1);
 DivisorX = CellDimension;
 DivisorY = CellDimension;
 DivisorZ = CellDimension;
 UCellLen = CellDimension*1000.0;
 
 fprintf(stderr,"Free Space Velocity     =%8.2E m/s\n",FreeSpaceVel );
 fprintf(stderr, "Operating frequency    = %g MHz \n", OperateFreq);
 fprintf(stderr, "Free Space Wave Number = %8.2f\n",WaveNumber);
 fprintf(stderr, "Free Space WaveLength  = %8.2f mm\n",WaveLength*1000.0);
 fprintf(stderr, "Unit cell dimension    = %10.4f  mm\n\n", UCellLen);
 fprintf(stderr, "Dimension of PEC box       = %g mm x %g mm x %g mm \n",
 Dimension_X*UCellLen,Dimension_Y*UCellLen,Dimension_Z*UCellLen);
 fprintf(stderr, "Total Number of cube cell  = %d  \n",TotNum_HexaHedron);
 fprintf(stderr, "Total Number of Nodes      = %d  \n",TotNum_GBL_Nodes);
}

/***************************************************************************
*
*   FUNCTION        : Assign_Global_Coordinates()
*   DESCRIPTION     : Assigns Global coordinates for each meshnode. 
*
****************************************************************************/

void Assign_Global_Coordinates()
 {
  int i,j,k,Xcount,Ycount,Zcount,Count_h=0; 
  double Xtemp,Ytemp,Ztemp;
   {
    Xcount=(int)Dimension_X;Ycount=(int)Dimension_Y;Zcount=(int)Dimension_Z;
    Xtemp=Min_X;  Ytemp=Min_Y;  Ztemp=Min_Z;
     for(k=0;k<=Zcount;k++) {
      for(j=0;j<=Ycount;j++) {
       for(i=0;i<=Xcount;i++) {
          MeshNodeCord[Count_h][0]=Xtemp; 
          MeshNodeCord[Count_h][1]=Ytemp; 
          MeshNodeCord[Count_h][2]=Ztemp; 
          Xtemp=Xtemp+DivisorX; Count_h++; }
          Xtemp=Min_X; Ytemp=Ytemp+DivisorY; }
          Ytemp=Min_Y; Ztemp=Ztemp+DivisorZ; }
    }
 }

/***************************************************************************
*
*   FUNCTION        : Assign_hexahedron_NodeNumbering()
*   DESCRIPTION     : Assigns Global nodenumbering for each meshnode.
*
****************************************************************************/

void Assign_hexahedron_NodeNumbering()
  {
   int i,start=1;
   {
   for(i=1;i<=TotNum_HexaHedron;i++)
     {
      HexHed[i-1][0]=start;
      HexHed[i-1][1]=HexHed[i-1][0]+1;
      HexHed[i-1][2]=HexHed[i-1][0]+Dimension_X+1;
      HexHed[i-1][3]=HexHed[i-1][2]+1;
      HexHed[i-1][4]=HexHed[i-1][0]+(Dimension_X+1)*(Dimension_Y+1);
      HexHed[i-1][5]=HexHed[i-1][4]+1;
      HexHed[i-1][6]=HexHed[i-1][4]+(Dimension_X+1);
      HexHed[i-1][7]=HexHed[i-1][6]+1;

      if((i%(Dimension_Y*Dimension_X)==0))
            { start=start+(Dimension_X+3); }
      else  { start=start+1; }

      if((start%(Dimension_X+1))==0)
           { start++; }

     }
  }
 }

/***************************************************************************
*
*   FUNCTION        : HexaHedra_SubDivision_1()
*   DESCRIPTION     : Divide each Hexahedron into five Tetrahedron.
*                     Find links between them.  (Method 1)
*
****************************************************************************/

void HexaHedra_SubDivision_1()
    {
      tetrnod[0][0]=HexHed[cn][0];   tetlink[0][0]=0;
      tetrnod[0][1]=HexHed[cn][2];   tetlink[0][1]=2;
      tetrnod[0][2]=HexHed[cn][4];   tetlink[0][2]=4;
      tetrnod[0][3]=HexHed[cn][1];   tetlink[0][3]=1;

      tetrnod[1][0]=HexHed[cn][5];   tetlink[1][0]=5;
      tetrnod[1][1]=HexHed[cn][1];   tetlink[1][1]=1;
      tetrnod[1][2]=HexHed[cn][4];   tetlink[1][2]=4;
      tetrnod[1][3]=HexHed[cn][7];   tetlink[1][3]=7;

      tetrnod[2][0]=HexHed[cn][3];   tetlink[2][0]=3;
      tetrnod[2][1]=HexHed[cn][1];   tetlink[2][1]=1;
      tetrnod[2][2]=HexHed[cn][7];   tetlink[2][2]=7;
      tetrnod[2][3]=HexHed[cn][2];   tetlink[2][3]=2;

      tetrnod[3][0]=HexHed[cn][6];   tetlink[3][0]=6;
      tetrnod[3][1]=HexHed[cn][2];   tetlink[3][1]=2;
      tetrnod[3][2]=HexHed[cn][7];   tetlink[3][2]=7;
      tetrnod[3][3]=HexHed[cn][4];   tetlink[3][3]=4;

      tetrnod[4][0]=HexHed[cn][1];   tetlink[4][0]=1;
      tetrnod[4][1]=HexHed[cn][2];   tetlink[4][1]=2;
      tetrnod[4][2]=HexHed[cn][4];   tetlink[4][2]=4;
      tetrnod[4][3]=HexHed[cn][7];   tetlink[4][3]=7;
   }

/***************************************************************************
*
*   FUNCTION        : HexaHedra_SubDivision_2()
*   DESCRIPTION     : Divide each Hexahedron into five Tetrahedron.
*                     Find links between them. (Method 2)
*
****************************************************************************/

void HexaHedra_SubDivision_2()
    {
      tetrnod[0][0]=HexHed[cn][2];   tetlink[0][0]=2;
      tetrnod[0][1]=HexHed[cn][0];   tetlink[0][1]=0;
      tetrnod[0][2]=HexHed[cn][3];   tetlink[0][2]=3;
      tetrnod[0][3]=HexHed[cn][6];   tetlink[0][3]=6;
 
      tetrnod[1][0]=HexHed[cn][7];   tetlink[1][0]=7;
      tetrnod[1][1]=HexHed[cn][6];   tetlink[1][1]=6;
      tetrnod[1][2]=HexHed[cn][3];   tetlink[1][2]=3;
      tetrnod[1][3]=HexHed[cn][5];   tetlink[1][3]=5;
 
      tetrnod[2][0]=HexHed[cn][1];   tetlink[2][0]=1;
      tetrnod[2][1]=HexHed[cn][0];   tetlink[2][1]=0;
      tetrnod[2][2]=HexHed[cn][5];   tetlink[2][2]=5;
      tetrnod[2][3]=HexHed[cn][3];   tetlink[2][3]=3;
 
      tetrnod[3][0]=HexHed[cn][4];   tetlink[3][0]=4;
      tetrnod[3][1]=HexHed[cn][5];   tetlink[3][1]=5;
      tetrnod[3][2]=HexHed[cn][0];   tetlink[3][2]=0;
      tetrnod[3][3]=HexHed[cn][6];   tetlink[3][3]=6;
 
      tetrnod[4][0]=HexHed[cn][3];   tetlink[4][0]=3;
      tetrnod[4][1]=HexHed[cn][0];   tetlink[4][1]=0;
      tetrnod[4][2]=HexHed[cn][5];   tetlink[4][2]=5;
      tetrnod[4][3]=HexHed[cn][6];   tetlink[4][3]=6;

   }

/***************************************************************************
*
*   FUNCTION        : Find_Volume_Determinant()
*   DESCRIPTION     : Find the Volume determinant of a Tetrahedron.
*
****************************************************************************/

void Find_Volume_Determinant()
     {
      int Count_i,Count_l,Count_j,Count_k,lu[4];
       {
       for(Count_l=0;Count_l<=3;Count_l++)
         {
         lu[Count_l]=tetrnod[t][Count_l];  
         }
       for(Count_i=0;Count_i<=3;Count_i++)
         for(Count_j=0;Count_j<=3;Count_j++)
           {
           if(Count_j==3)
            {
            Volume_Cord[Count_i][Count_j]=1.0;
            VVolume_Cord[Count_i][Count_j]=1.0;
            }
           else
            {
            Count_k=lu[Count_i];
            Volume_Cord[Count_i][Count_j]=MeshNodeCord[Count_k-1][Count_j];
            VVolume_Cord[Count_i][Count_j]=MeshNodeCord[Count_k-1][Count_j];
            } 
           }
         {
         det=Determinant_TetraMatrix();
         }
      }
    }

/***************************************************************************
*
*   FUNCTION        : Find_CoFactor() 
*   DESCRIPTION     : Find the  Cofactors of a Volume determinant.
*
****************************************************************************/

void Find_CoFactor()
   {
   short int i,j,row,col,ro,co,counter=0;
   double cof[4][4];
     {

     for(row=0;row<=3;row++)
     for(col=0;col<=3;col++) {
     co=0;ro=0;counter=0;  

        for(i=0;i<=3;i++)
        for(j=0;j<=3;j++) {
        if((i!=row) && (j!=col)) {
        counter++;
        if(counter<4) ro=0;
        if((counter>3)&&(counter<7)) ro=1;
        else if(counter>6) ro=2;
        cof[ro][co]=VVolume_Cord[i][j];
        co++;co=co%3;            }
                          }
        ro=0;co=0;counter=0;
        cofs[row][col]=cof[0][0]*((cof[1][1]*cof[2][2])-(cof[1][2]*cof[2][1]))
                      -cof[0][1]*((cof[1][0]*cof[2][2])-(cof[2][0]*cof[1][2]))
                      +cof[0][2]*((cof[1][0]*cof[2][1])-(cof[2][0]*cof[1][1]));
        
        if((row+col)%2!=0)
        cofs[row][col]*=-1.0;
         }
      }
    }

/**************************************************************************
*
*   FUNCTION        : TetraHedron_SubMatrix()
*   DESCRIPTION     : Form 12 x 12 Submatrix for a Tetrahedron.
*
**************************************************************************/

void TetraHedron_SubMatrix()
{
 short int mm0,mm1,mm2,mm3,cnt1,cnt2,p,q,m,n,i,k;
 double real_perm,real_buff,delta;
 { 

   mm0=tetrnod[t][0];
   mm1=tetrnod[t][1];
   mm2=tetrnod[t][2];
   mm3=tetrnod[t][3];
   if((REL_Permittivity[mm0]==1.0) && (REL_Permittivity[mm1]==1.0) &&
      (REL_Permittivity[mm2]==1.0) && (REL_Permittivity[mm3]==1.0))
    {
        for(m=1;m<=4;m++)
        for(n=1;n<=4;n++)
        for(i=1;i<=3;i++)                          {
        p=3*(m-1)+i; q=3*(n-1)+i;
        if(m==n) delta=1.0; else delta=0.0;
        real_buff=((1+delta)*det)/20.0;
        real_buff=real_buff*(WaveNumber*WaveNumber);
        Tet_Matrix_Element[t][p-1][q-1]=real_buff; }
    }
   else
    {
        for(m=1;m<=4;m++)
        for(n=1;n<=4;n++)
          {
           cnt1=tetrnod[t][m-1]; cnt2=tetrnod[t][n-1];
           if(REL_Permittivity[cnt1]==REL_Permittivity[cnt2])  
             {
                real_perm=REL_Permittivity[cnt1];
                for(i=1;i<=3;i++)                          {
                p=3*(m-1)+i; q=3*(n-1)+i;
                if(m==n) delta=1.0; else delta=0.0;
                real_buff=((1+delta)*det)/20.0;
                real_buff=real_buff*(WaveNumber*WaveNumber);
                if(REL_Permittivity[cnt1]==1.0)
                     { if(i<=2) {real_buff=real_buff*real_perm;}}
                else {           real_buff=real_buff*real_perm; }
                Tet_Matrix_Element[t][p-1][q-1]=real_buff; }
             }
           else
             {
                if(REL_Permittivity[cnt1]>=REL_Permittivity[cnt2]) 
                  {real_perm=REL_Permittivity[cnt1];}
                else
                  {real_perm=REL_Permittivity[cnt2];}

                for(i=1;i<=3;i++)                           {
                p=3*(m-1)+i; q=3*(n-1)+i;
                if(m==n) delta=1.0; else delta=0.0;
                real_buff=((1+delta)*det)/20.0;
                real_buff=real_buff*(WaveNumber*WaveNumber);
                if(i<=2) {real_buff=real_buff*real_perm;}
                Tet_Matrix_Element[t][p-1][q-1]=real_buff;  }
              }

          } 

     }
  for(m=1;m<=4;m++)          
    for(n=1;n<=4;n++)
      for(i=1;i<=3;i++)                          
        {
        p=3*(m-1)+i; q=3*(n-1)+i;
        real_buff=(cofs[m-1][0]*cofs[n-1][0])+
                  (cofs[m-1][1]*cofs[n-1][1])+
                  (cofs[m-1][2]*cofs[n-1][2]);
        real_buff=(real_buff)/(36.0*det);
        Tet_Matrix_Element[t][p-1][q-1]-=real_buff;
        }

  for(m=1;m<=4;m++)
    for(n=1;n<=4;n++)
      for(i=1;i<=3;i++)
        for(k=1;k<=3;k++)
        {
        if(i!=k) {
          p=3*(m-1)+i; q=3*(n-1)+k;
          real_buff=(cofs[m-1][k-1]*cofs[n-1][i-1]);
          real_buff-=(cofs[m-1][i-1]*cofs[n-1][k-1]);
          real_buff=(real_buff)/(36.0*det);
          Tet_Matrix_Element[t][p-1][q-1]+=real_buff; }
        }

 for(p=0;p<=11;p++) for(q=0;q<=11;q++)
 if(fabs(Tet_Matrix_Element[t][p][q]) < 1.0e-10) 
 {Tet_Matrix_Element[t][p][q] = 0.0;}

 }
}

/***************************************************************************
*
*   FUNCTION        : HexaHedron_Matrix()
*   DESCRIPTION     : Form 24 x 24 Submatrix for a Hexahedron.
*
****************************************************************************/

void HexaHedron_Matrix()
   {
    short int i,j,t,mm,nn,m,n,k,l;
    {

   for(i=0;i<=23;i++) for(j=0;j<=23;j++)
   { Hex_Matrix_Element[i][j]=0.0;}
 
   for(t=0;t<=4;t++) 
    for(m=0;m<=3;m++) 
      for(n=0;n<=3;n++)
        {
          k=tetlink[t][m]*3; 
          l=tetlink[t][n]*3;
          mm=m*3; nn=n*3;
          for(i=0;i<=2;i++) 
            for(j=0;j<=2;j++)
              {
               if(Tet_Matrix_Element[t][mm+i][nn+j]!=0.0)
               {
               Hex_Matrix_Element[k+i][l+j]+=
               Tet_Matrix_Element[t][mm+i][nn+j];
               }
              }
        }
  for(m=0;m<=23;m++) for(n=0;n<=23;n++)
  if(fabs(Hex_Matrix_Element[m][n]) < 1.0e-10)
   {Hex_Matrix_Element[m][n] = 0.0;}
    }
   }

/***************************************************************************
*
*   FUNCTION        : HexaHedron_link_1()
*   DESCRIPTION     : Provide links between the nodes. ( Method 1)
*
****************************************************************************/

void HexaHedron_link_1()
    {
     connect1[0][0] = 1;  connect1[0][1] = 2;
     connect1[0][2] = 4;  connect1[0][3] =-1;
     connect1[0][4] =-1;  connect1[0][5] =-1;

     connect1[1][0] = 0;  connect1[1][1] = 2;
     connect1[1][2] = 3;  connect1[1][3] = 4;
     connect1[1][4] = 5;  connect1[1][5] = 7;

     connect1[2][0] = 0;  connect1[2][1] = 1;
     connect1[2][2] = 3;  connect1[2][3] = 4;
     connect1[2][4] = 6;  connect1[2][5] = 7;

     connect1[3][0] = 1;  connect1[3][1] = 2;
     connect1[3][2] = 7;  connect1[3][3] =-1;
     connect1[3][4] =-1;  connect1[3][5] =-1;

     connect1[4][0] = 0;  connect1[4][1] = 1;
     connect1[4][2] = 2;  connect1[4][3] = 5;
     connect1[4][4] = 6;  connect1[4][5] = 7;

     connect1[5][0] = 1;  connect1[5][1] = 4;
     connect1[5][2] = 7;  connect1[5][3] =-1;
     connect1[5][4] =-1;  connect1[5][5] =-1;

     connect1[6][0] = 2;  connect1[6][1] = 4;
     connect1[6][2] = 7;  connect1[6][3] =-1;
     connect1[6][4] =-1;  connect1[6][5] =-1;

     connect1[7][0] = 1;  connect1[7][1] = 2;
     connect1[7][2] = 3;  connect1[7][3] = 4;
     connect1[7][4] = 5;  connect1[7][5] = 6;
     }

/***************************************************************************
*
*   FUNCTION        : HexaHedron_link_2()
*   DESCRIPTION     : Provide links between the nodes. ( Method 2)
*
****************************************************************************/
 
void HexaHedron_link_2()
    {
     connect1[0][0] = 1;  connect1[0][1] = 2;
     connect1[0][2] = 3;  connect1[0][3] = 4;
     connect1[0][4] = 5;  connect1[0][5] = 6;
 
     connect1[1][0] = 0;  connect1[1][1] = 3;
     connect1[1][2] = 5;  connect1[1][3] =-1;
     connect1[1][4] =-1;  connect1[1][5] =-1;
 
     connect1[2][0] = 0;  connect1[2][1] = 3;
     connect1[2][2] = 6;  connect1[2][3] =-1;
     connect1[2][4] =-1;  connect1[2][5] =-1;
 
     connect1[3][0] = 0;  connect1[3][1] = 1;
     connect1[3][2] = 2;  connect1[3][3] = 5;
     connect1[3][4] = 6;  connect1[3][5] = 7;
 
     connect1[4][0] = 0;  connect1[4][1] = 5;
     connect1[4][2] = 6;  connect1[4][3] =-1;
     connect1[4][4] =-1;  connect1[4][5] =-1;
 
     connect1[5][0] = 0;  connect1[5][1] = 1;
     connect1[5][2] = 3;  connect1[5][3] = 4;
     connect1[5][4] = 6;  connect1[5][5] = 7;
 
     connect1[6][0] = 0;  connect1[6][1] = 2;
     connect1[6][2] = 3;  connect1[6][3] = 4;
     connect1[6][4] = 5;  connect1[6][5] = 7;
 
     connect1[7][0] = 3;  connect1[7][1] = 5;
     connect1[7][2] = 6;  connect1[7][3] =-1;
     connect1[7][4] =-1;  connect1[7][5] =-1;
     }

/***************************************************************************
*
*   FUNCTION        : Find_Global_Matrix()
*   DESCRIPTION     : Construct 3N x 3N Global matrix. N =Total # of nodes.
*
****************************************************************************/

void Find_Global_Matrix()
{
 int l,s,x,y,z,i,j,ma,qa,k=0,c;
 int tem=0,rval,fval,bck;
{

for(k=0;k<=7;k++)
{

ma=HexHed[cn][k]; ma--;s=ma*3;c=k*3; 

for(i=0;i<=2;i++)
 for(j=0;j<=2;j++)   
  {
  if(fabs(Hex_Matrix_Element[c+i][c+j])>=1.0e-10)
  {
   tem=GBL_Matrix_Index[s+i];
   {rval=Srch_NZ_Element_Column(s+i,s+j);}
   if(rval>=0)
   {
    bck=tem;tem=rval;
    GBL_Matrix_ColNos[s+i][tem]=s+j;
    GBL_Matrix_Data[s+i][tem]+=Hex_Matrix_Element[c+i][c+j];
    tem=bck;
    if(tem==0) tem++; 
    GBL_Matrix_Index[s+i]=tem;
   }
   else
   {
   GBL_Matrix_ColNos[s+i][tem]=s+j;
   GBL_Matrix_Data[s+i][tem]=Hex_Matrix_Element[c+i][c+j];
   tem++;
   GBL_Matrix_Index[s+i]=tem;
   }
  }
  }

for(y=0;y<=5;y++)
{
 x=connect1[k][y];  if(x==-1) continue; 

 l=HexHed[cn][x]; l--; z=l*3;qa=x*3;

 for(i=0;i<=2;i++)
  for(j=0;j<=2;j++)
   {
   if(Hex_Matrix_Element[c+i][qa+j]!=0.0)
   {
   tem=GBL_Matrix_Index[s+i];
   {fval=Srch_NZ_Element_Column(s+i,z+j);}
   if(fval>=0)
   {
    bck=tem;tem=fval;
    GBL_Matrix_ColNos[s+i][tem]=z+j;
    GBL_Matrix_Data[s+i][tem]+=Hex_Matrix_Element[c+i][qa+j];
    tem=bck;
    GBL_Matrix_Index[s+i]=tem;
   }
   else
   {
   GBL_Matrix_ColNos[s+i][tem]=z+j;
   GBL_Matrix_Data[s+i][tem]=Hex_Matrix_Element[c+i][qa+j];
   tem++;
   GBL_Matrix_Index[s+i]=tem;
   }
  }
 }
}
ma=0;qa=0;c=0;s=0;z=0;x=0;
}
}
}

/***************************************************************************
*
*   FUNCTION        : Read_Input_Pass_1()
*   DESCRIPTION     : find the number of nodes and the number of hexahedron.
*
****************************************************************************/
 
void Read_Input_Pass_1()
{
     int    x1, y1, z1, x2, y2, z2, OutF_Count=0, output_flag=0;
     char   type[20], att1[20], att2[18], att3[18], att4[18];
     double tmp1, tmp2, cell_dim=1.0;
     char   buffer[80];
 
     /* READ THE INPUT FILE:  FIRST PASS */
 
        while (fgets(buffer, 80, InF)) {
        if (buffer[0] == '#') {continue;}
        if (!strncmp(buffer, "celldim", 7)) { sscanf(buffer, "%s%lf%s", type,
                                              &tmp1, att1);
                                              if (!strcmp(att1, "cm")) tmp2=0.01;
                                              else if(!strcmp(att1, "m"))
                                                      tmp2=1.0;
                                              else if(!strcmp(att1, "mm"))
                                                      tmp2=0.001;
                                              else if(!strcmp(att1, "in"))
                                                      tmp2=2.54*0.01;
                                              else {fprintf(stderr, "unrecognized cell dimension.\n");
                                                    exit(1);}
                                              cell_dim=tmp1*tmp2;
                                              continue;
                                             }
        if(!strncmp(buffer, "default_output",14))
                   {
                   sscanf(buffer, "%s%s", type, Out_FileName0);
                   OutF_0=fopen(Out_FileName0, "w");
                   if (OutF_0 == NULL) {fprintf(stderr, " Error: Can't create default output file %s\n", Out_FileName0); exit(1);}
                   output_flag=1;
                   continue;
                   }
        else if(sscanf(buffer,"%s%d%d%d%d%d%d%s%s%s%s",type,&x1,&y1,&z1,&x2,&y2,&z2,att1,att2,att3,att4)>2)
           {
             if(!strcmp(type, "boundary"))
             {
/*
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   if((x1==x2)||(y1==y2)||(z1==z2)) {fprintf(stderr,"ERROR: boundary must have finite thickness/n");
                                                     exit(1);
                                                    }
                   continue;
*/
fprintf(stderr, "\n\nERROR: boundary must be specified by the option 'box' \n"); exit(1);
             }
             if(!strcmp(type, "dielectric"))
             {
                   /* fprintf(stderr, "dielectric \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   if((x1==x2)||(y1==y2)||(z1==z2)){fprintf(stderr,"ERROR: dielectric must have finite thickness/n");
                                                     exit(1);
                                                    }
                   continue;
             }
             if(!strcmp(type, "box"))
             {
                   /* fprintf(stderr, "box \n"); */
                   if((x1==x2)||(y1==y2)||(z1==z2)){ fprintf(stderr, "Error: invalid box:  ");
                                                   fprintf(stderr, "box must be three dimensional\n");
                                                   exit(1);}
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   continue;
             }
             if(!strcmp(type,"conductor"))
             {
                   /* fprintf(stderr, "conductor \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   continue;
             }
             if(!strcmp(type, "aperture"))
             {
                   /* fprintf(stderr, "aperture \n"); */
                if((x1!=x2) & (y1!=y2) & (z1!=z2)) {fprintf(stderr, "ERROR: apertures must be two dimensional\n");
                                                       exit(1);
                                                      }
                if(((x1==x2)&&(y1==y2))||((x1==x2)&&(z1==z2))||((z1==z2)&&(y1==y2)))
                                                   {fprintf(stderr, "ERROR: apertures must be two dimensional\n");
                                                       exit(1);
                                                      }
                   continue;
             }
             if(!strcmp(type, "seam"))
             {
                   /* fprintf(stderr, "seam \n"); */
                   continue;
             }
             if(!strcmp(type,"esource"))
             {
                   /* fprintf(stderr, "esource \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   freq = atof(att1);
                   if((x1==x2)&&(y1==y2)&&(z1==z2)) {
                                  fprintf(stderr, "ERROR: esource cannot be defined at a point\n");
                                  exit(1);
                                  }
                   if ((att2[0] != 'x')  && (att2[0] != 'y') && (att2[0] != 'z')){
                         fprintf(stderr, "ERROR: unrecognized polarization of esource statement  %c \n", att2[0]);
                                  exit(1);
                                  }
                   continue;
              }
              if(!strcmp(type, "msource"))
              {
                   /* fprintf(stderr, "msource \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   freq = atof(att1);
                   continue;
              }
              if(!strcmp(type, "efield_output"))
              {    
                   switch (OutF_Count)
                        {
                        case 0:  {OutF_Count=1; strcpy(Out_FileName1, att1); OutF_1=fopen(Out_FileName1,"w"); break;}
                        case 1:  {if(strcmp(att1, Out_FileName1)){OutF_Count=2; strcpy(Out_FileName2, att1); OutF_2=fopen(Out_FileName2,"w");}
                                  break;
                                 }
                        case 2:  {if(strcmp(att1, Out_FileName1) && strcmp(att1, Out_FileName2))
                                      {OutF_Count=3; strcpy(Out_FileName3, att1); OutF_3=fopen(Out_FileName3,"w");}
                                  break;
                                 }
                        case 3:  {if((strcmp(att1, Out_FileName1) && strcmp(att1, Out_FileName2)) && strcmp(att1, Out_FileName3))
                                      {fprintf(stderr, "ERROR: too many different files specified by efield_output commands \n"); exit(1);}
                                  break;
                                 }
                        }
                   output_flag=1;
                   continue;
              }
              fprintf(stderr, "Unrecognized attribute %s\n\n", type);
            }                 /*end of else*/
        }                     /* end of while */
 
 
        if(freq ==0) fprintf(stderr, "\nWARNING: No sources were specified, code will not produce output!\n\n");
        if(output_flag == 0) fprintf(stderr,"\nWARNING: The input file does not contain an output statement, code will not produce output!\n\n");
        OperateFreq = freq;  CellDimension = cell_dim;
        Dimension_X=Max_X-Min_X; Dimension_Y=Max_Y-Min_Y; Dimension_Z=Max_Z-Min_Z;
        TotNum_HexaHedron=Dimension_X*Dimension_Y*Dimension_Z;
        TotNum_GBL_Nodes=(Dimension_X+1)*(Dimension_Y+1)*(Dimension_Z+1);
        rewind(InF);
 }

/***************************************************************************
*
*   FUNCTION        : Read_Input_Pass_2()
*   DESCRIPTION     : Reads information for each node from input datafile.
*                     Reads information about boundary conditions.
*
****************************************************************************/

void Read_Input_Pass_2()
  {
     int    i, j, k, I, J, x1, y1, z1, x2, y2, z2, nn, tmp;
     char   type[20], att1[20], att2[18], att3[18], att4[18];
     char   buffer[80];

  /*  Initialize dynamically allocated variables  */
    for(nn=0;nn<TotNum_GBL_Nodes;nn++)
       {for(i=0;i<=2;i++){Dat_Buffer[nn][i]=-999.;REL_Permittivity[nn]=1.0;}}

     I=Dimension_X+1; J=Dimension_Y+1;

   /*Begin reading the input file for the second time */
        while (fgets(buffer, 80, InF))
        {
            if (buffer[0] == '#') continue;
            if(sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s", 
            type, &x1, &y1, &z1, &x2, &y2, &z2, att1, att2, att3, att4)>2)
             {
             if(!strcmp(type, "boundary"))
             {
                   continue;
             }
             if(!strcmp(type, "efield_output"))
             {     
                   if((x1<Min_X || x1>Max_X)||(x2<Min_X || x2>Max_X)) {fprintf(stderr,"ERROR: efield_output parameter out of range.\n"); exit(1);}
                   if((y1<Min_Y || y1>Max_Y)||(y2<Min_Y || y2>Max_Y)) {fprintf(stderr,"ERROR: efield_output parameter out of range.\n"); exit(1);}
                   if((z1<Min_Z || z1>Max_Z)||(z2<Min_Z || z2>Max_Z)) {fprintf(stderr,"ERROR: efield_output parameter out of range.\n"); exit(1);}
                   continue;
             }

             if(!strcmp(type, "box"))
             {
                   if (x2<x1) swap(&x1, &x2); if (y2<y1) swap(&y1, &y2); if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;

                   for(j=y1;j<=y2;j++){
                     for(k=z1;k<=z2;k++){
                                         nn=k*J*I+j*I+x1;
                                         Dat_Buffer[nn][1]=0.0;
                                         Dat_Buffer[nn][2]=0.0;
                                         nn=k*J*I+j*I+x2;
                                         Dat_Buffer[nn][1]=0.0;
                                         Dat_Buffer[nn][2]=0.0;
                                        }
                                      }                           
                   for(i=x1;i<=x2;i++){
                     for(k=z1;k<=z2;k++){
                                         nn=k*J*I+y1*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][2]=0.0;
                                         nn=k*J*I+y2*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][2]=0.0;                              
                                        }
                                      }                           
                   for(i=x1;i<=x2;i++){
                     for(j=y1;j<=y2;j++){
                                         nn=z1*J*I+j*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][1]=0.0;
                                         nn=z2*J*I+j*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][1]=0.0;
                                        }
                                      }                           
                  continue;
             }
             if(!strcmp(type, "dielectric"))
             {
                   if (x2<x1) swap(&x1, &x2); if (y2<y1) swap(&y1, &y2); if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;
                   for(k=z1;k<=z2;k++){
                      for(j=y1;j<=y2;j++){
                         for(i=x1;i<=x2;i++){
                         nn=k*I*J+j*I+i;
                         REL_Permittivity[nn]=atof(att1);
                                            }
                                         }
                                      } 
                   continue;
             }
             if(!strcmp(type, "conductor"))
             {
                   if (x2<x1) swap(&x1, &x2); if (y2<y1) swap(&y1, &y2); if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;

              /* three-dimensional conductor */
                   for(i=x1;i<=x2;i++){
                     for(j=y1;j<=y2;j++){
                        for(k=z1;k<=z2;k++){
                                           nn=k*J*I+j*I+i;
                                           for(tmp=0;tmp<2;tmp++){
                                                                  Dat_Buffer[nn][tmp]=0.0;
                                                                 }
                                          }
                                       }
                                     }
               /* two-dimensional conductor */
                  if(x1==x2){
                           for(j=y1;j<=y2;j++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+j*I+x1;
                                                Dat_Buffer[nn][1]=0.0;
                                                Dat_Buffer[nn][2]=0.0;
                                               }
                                              }                           
                           }
                 if(y1==y2){
                           for(i=x1;i<=x2;i++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+y1*I+i;
                                                Dat_Buffer[nn][0]=0.0;
                                                Dat_Buffer[nn][2]=0.0;                                                
                                               }
                                              }                           
                           }
                 if(z1==z2){
                           for(i=x1;i<=x2;i++){
                             for(j=y1;j<=y2;j++){
                                                nn=z1*J*I+j*I+i;
                                                Dat_Buffer[nn][0]=0.0;
                                                Dat_Buffer[nn][1]=0.0;
                                               }
                                              }                           
                           }
               /* one-dimensional conductor */
                  if ((x1==x2)&&(y1==y2)){
                                          for(k=z1;k<z2;k++){
                                                             nn=k*J*I+y1*I+x1;
                                                             Dat_Buffer[nn][2]=0.0;
                                                            } 
                                         } 
                  if ((x1==x2)&&(z1==z2)){
                                          for(j=y1;j<y2;j++){
                                                             nn=z1*J*I+j*I+x1;
                                                             Dat_Buffer[nn][1]=0.0;
                                                            } 
                                         }
                  if ((y1==y2)&&(z1==z2)){
                                          for(i=x1;i<x2;i++){
                                                             nn=z1*J*I+y1*I+i;
                                                             Dat_Buffer[nn][0]=0.0;
                                                            } 
                                         }                                           
                   continue;
             }
             if(!strcmp(type, "aperture"))
             {
                   /* fprintf(stderr, "aperture \n"); */
                  if(x1==x2){
                           for(j=y1;j<=y2;j++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+j*I+x1;
                                                Dat_Buffer[nn][0]=-999.;
                                                Dat_Buffer[nn][1]=-999.;
                                                Dat_Buffer[nn][2]=-999.;
                                               }
                                              }                           
                           }
                 if(y1==y2){
                           for(i=x1;i<=x2;i++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+y1*I+i;
                                                Dat_Buffer[nn][0]=-999.;
                                                Dat_Buffer[nn][1]=-999.;
                                                Dat_Buffer[nn][2]=-999.;                                                
                                               }
                                              }                           
                           }
                 if(z1==z2){
                           for(i=x1;i<=x2;i++){
                             for(j=y1;j<=y2;j++){
                                                nn=z1*J*I+j*I+i;
                                                Dat_Buffer[nn][0]=-999.;
                                                Dat_Buffer[nn][1]=-999.;
                                                Dat_Buffer[nn][2]=-999.;
                                               }
                                              }                           
                           }

                   continue;
             }
             if(!strcmp(type, "seam"))
             {
                   /* fprintf(stderr, "seam \n"); */
                   /*  need to define seam here  */
                   continue;
             }
             if(!strcmp(type, "esource"))
             {
                   /* fprintf(stderr, "esource \n"); */
                   if(freq != atof(att1)) { fprintf(stderr, "Error: all sources must have the same frequency\n");
                                            exit(1);
                                          }
                   if (x2<x1) swap(&x1, &x2); if (y2<y1) swap(&y1, &y2); if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;

                  /* three dimensional source blocks */       
                   if (att2[0] == 'y')
                     {
                      for(i=x1;i<=x2;i++){
                       for(j=y1;j<=y2;j++){
                        for(k=z1;k<=z2;k++){
                                            nn =k*I*J+j*I+i;
                                            Dat_Buffer[nn][1]=atof(att3);
                                           }
                                          }
                                         }
                     }

                   if (att2[0] == 'z')
                     { 
                      for(i=x1;i<=x2;i++){
                       for(j=y1;j<=y2;j++){
                        for(k=z1;k<=z2;k++){
                                           nn = k*I*J+j*I+i;
                                           Dat_Buffer[nn][2]=atof(att3);
                                          }
                                         }
                                        }
                     }

                  if (att2[0] == 'x')
                     {
                      for(i=x1;i<=x2;i++){
                       for(j=y1;j<=y2;j++){
                        for(k=z1;k<=z2;k++){
                                           nn = k*I*J+j*I+i;
                                           Dat_Buffer[nn][0]=atof(att3);
                                          }
                                         }
                                        }
                     }

                 continue;
              }
              if(!strcmp(type, "ignore")) continue;
              fprintf(stderr, "Warning: Unrecognized attribute - %s\n", type);
           }              /*end of if*/ 
        }                 /* end of while */
       rewind(InF);
       /* Count the number of known and unknown vectors */
       TotNum_Known_Var=0;
       TotNum_Unknown_Var=0;
       for(nn=0;nn<TotNum_GBL_Nodes;nn++){
           for(i=0;i<=2;i++){
            if(Dat_Buffer[nn][i]==-999.) TotNum_Unknown_Var++;
            else                         TotNum_Known_Var++;
                            }                                                                                          
                                         } 
  }

/***************************************************************************
*
*   FUNCTION        : Fill_Vectors()
*   DESCRIPTION     : Defines forced and unforced vectors
*
****************************************************************************/

void Fill_Vectors()

  {
       int    i, nn;
 
       /* Define forced and unforced vectors */
       TotNum_Known_Var=0;
       TotNum_Unknown_Var=0;
       for(nn=0;nn<TotNum_GBL_Nodes;nn++){
           for(i=0;i<=2;i++){
                             if(Dat_Buffer[nn][i]==-999.){
                                                         UN_Vector[TotNum_Unknown_Var]=nn*3+i;
                                                         Loc_UN_Node[TotNum_Unknown_Var]=nn;
                                                         Loc_UN_NodeCord[TotNum_Unknown_Var]=i;
                                                         TotNum_Unknown_Var++;
                                                        }
                             else if(Dat_Buffer[nn][i]!=0.0){
                                                             FORC_Vector[TotNum_Known_Var]=nn*3+i;
                                                             FORC_Val[TotNum_Known_Var]=Dat_Buffer[nn][i];
                                                             TotNum_Known_Var++;
                                                            }                                                                                          
                            }
                                         } 
  }

/***************************************************************************
*
*   FUNCTION        : Count_Half_BandWidth()
*   DESCRIPTION     : Counts Bandwidth of the matrix A of Ax = B.
*
****************************************************************************/

 void Count_Half_BandWidth()
     {
      short int i,j,k,l,ColDiff,Val;
        {
        for(i=0;i<=TotNum_Unknown_Var-1;i++)
          {
          k=UN_Vector[i];
          for(j=i;j<=TotNum_Unknown_Var-1;j++)
            {
            l=UN_Vector[j];
            Val=Srch_NZ_Element_Column(k,l);
            if((Val>=0))
             {
              ColDiff=j-i;
              if(Half_BandWidth<=(ColDiff+1))
               {
                Half_BandWidth=(ColDiff+1);
               }
             }
           }
         }
         printf("Half_BandWidth = %d \n",Half_BandWidth);
      }
    }

/***************************************************************************
*
*   FUNCTION      : Partition_Global_Matrix()
*   DESCRIPTION   : Partitions Global matrix and applies boundary conditions.
*                   Construct A matrix of Ax = B matrix equation.
*
****************************************************************************/

void Partition_Global_Matrix()
     {
      short int i,j,k,l,ColDiff,g=0,val;
        {

        for(i=0;i<=TotNum_Unknown_Var-1;i++)  { k=UN_Vector[i];
        for(j=i;j<=TotNum_Unknown_Var-1;j++)  { l=UN_Vector[j];
        val=Srch_NZ_Element_Column(k,l);
        if((val>=0))
         {
          ColDiff=j-i;
          HalfBand_Matrix[i+1][ColDiff+1]+=GBL_Matrix_Data[k][val];
         }
                                              }

         g=0;
         for(j=0;j<=TotNum_Known_Var-1;j++)  { l=FORC_Vector[j];
         val=Srch_NZ_Element_Column(k,l);
         if(val>=0)
          {
           FORC_Matrix_ColNos[i][g]=j;
           FORC_Matrix_Data[i][g]=(GBL_Matrix_Data[k][val]);
           g++;
           }
                                             }
          g=0;                               }
          g=0;
      }
    }

/***************************************************************************
*
*   FUNCTION        : Find_RHS_Vector()
*   DESCRIPTION     : Calculates RHS vectors B of Ax = B matrix equation.
*
****************************************************************************/

void Find_RHS_Vector()
    {
     int Count_i,Count_j,Srch_Ptr=0; double Temp_Buff,Multpl;
     {
      for(Count_i=0;Count_i<=TotNum_Unknown_Var-1;Count_i++)
        {
        for(Count_j=0;Count_j<=TotNum_Known_Var-1;Count_j++)
          {
           Srch_Ptr=Locate_NZ_Element_Column(Count_i,Count_j);
           if(Srch_Ptr>=0) { Temp_Buff=FORC_Matrix_Data[Count_i][Srch_Ptr];}
           else            { Temp_Buff=0.0;} 
           if(Temp_Buff!=0.0)
            {
             Multpl=Temp_Buff*FORC_Val[Count_j];
             RHS_Vector[Count_i]+=Multpl;
            }
          }
          RHS_Vector[Count_i]=RHS_Vector[Count_i]*-1.0;
        }
      }
    }

/***************************************************************************
*
*   FUNCTION        : Produce_Output()
*   DESCRIPTION     : Prints the efield data to the output files.
*
****************************************************************************/

void Produce_Output()
{
 int       i, j, k, x1, y1, z1, x2, y2, z2, Count_i, dflag=0;
 short int ColNos, RowNos;
 char      type[20], att1[20], att2[8], att3[8], att4[8];
 char      buffer[80];
 FILE      *OutF;

  for(Count_i=0;Count_i<=TotNum_Unknown_Var-1;Count_i++)   /*put data in Dat_Buffer*/
   {
    RowNos=Loc_UN_Node[Count_i]; ColNos=Loc_UN_NodeCord[Count_i];
    Dat_Buffer[RowNos][ColNos]=RHS_Vector[Count_i];
   }

        while (fgets(buffer, 80, InF))
        {
         if (buffer[0] == '#') {
                               if(OutF_0 != NULL) fprintf(OutF_0, "%s", buffer);
                               if(OutF_1 != NULL) fprintf(OutF_1, "%s", buffer);
                               if(OutF_2 != NULL) fprintf(OutF_2, "%s", buffer);
                               if(OutF_3 != NULL) fprintf(OutF_3, "%s", buffer);
                               continue;
                              }
        if(!strncmp(buffer, "default_output",14)) {dflag=1; continue;}
        if(!strncmp(buffer, "efield_output",13))
                   { 
                   sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s", type, &x1, &y1, &z1, &x2, &y2, &z2, att1, att2, att3, att4);
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   if (!strcmp(att1, Out_FileName1)) OutF=OutF_1;
                   if (!strcmp(att1, Out_FileName2)) OutF=OutF_2;
                   if (!strcmp(att1, Out_FileName3)) OutF=OutF_3;
                   for(k=z1;k<=z2;k++){ 
                      for(j=y1;j<=y2;j++){
                         for(i=x1;i<=x2;i++){
                                             EfieldatNode(i, j, k, OutF); 
                                            }
                                         }
                                      } 
                   continue;
                   }

        }                 /* end of while */
        if(dflag == 1)    /* prints default output data */
              {

                for(Count_i=0;Count_i<=TotNum_GBL_Nodes-1;Count_i++)
                 {
                 fprintf(OutF_0," %6.3f\t%6.3f\t%6.3f\t",MeshNodeCord[Count_i][0],
                 MeshNodeCord[Count_i][1],MeshNodeCord[Count_i][2]);
                 fprintf(OutF_0," %f\t%f\t%f\t",
                 Dat_Buffer[Count_i][0],Dat_Buffer[Count_i][1],Dat_Buffer[Count_i][2]);
                 fprintf(OutF_0,"\n");
                 }
              }
        fclose(OutF_1);
        fclose(OutF_2);
        fclose(OutF_3);
        fclose(InF);
        fprintf(stdout,"\nProgram Successfully Completed: \n\n");
}

/***************************************************************************
*
*   FUNCTION        : Band_Matrix_Solver()
*   DESCRIPTION     : Solves for x variables of Ax = B matrix equation.
*
****************************************************************************/

void Band_Matrix_Solver()
  {
  int Count_i,Count_j,Count_k,Count_l,Count_m,Count_n;
  int Row_Number,Column_Number;
  double Pivot_Emement,Multpl_1,Multpl_2;
   {

 Row_Number=TotNum_Unknown_Var;
 Column_Number=Half_BandWidth;
 for(Count_i=1;Count_i<=Row_Number;Count_i++)
  {
  for(Count_j=2;Count_j<=Column_Number;Count_j++)
   {
   if(HalfBand_Matrix[Count_i][Count_j]!=0.0)
   {
   Pivot_Emement=HalfBand_Matrix[Count_i][Count_j]/HalfBand_Matrix[Count_i][1];
   Count_n=Count_i+Count_j-1;
   Count_l=0;
   for(Count_k=Count_j;Count_k<=Column_Number;Count_k++)
    {
    Count_l=Count_l+1;
    Multpl_1=Pivot_Emement*HalfBand_Matrix[Count_i][Count_k];
    HalfBand_Matrix[Count_n][Count_l] -= Multpl_1;
    }
    Multpl_2=Pivot_Emement*RHS_Vector[Count_i-1];
    RHS_Vector[Count_n-1] -= Multpl_2;
   }
  }
  }


 RHS_Vector[Row_Number-1] /= HalfBand_Matrix[Row_Number][1];
 for(Count_m=2;Count_m<=Row_Number;Count_m++)
  {
  Count_i=Row_Number+1-Count_m;
  for(Count_j=2;Count_j<=Column_Number;Count_j++)
   {
   if(HalfBand_Matrix[Count_i][Count_j]!=0.0)
   {
   Count_k=Count_i+Count_j-1;
   Multpl_1=HalfBand_Matrix[Count_i][Count_j]*RHS_Vector[Count_k-1];
   RHS_Vector[Count_i-1] -= Multpl_1;
   }
  }
  RHS_Vector[Count_i-1] /= HalfBand_Matrix[Count_i][1];
  }
 }
 }

/***************************************************************************
*
*   FUNCTION    : Srch_NZ_Element_Column()
*   DESCRIPTION : Searches for column number of non-zero element
*                 within each row of Global matrix. 
*
****************************************************************************/

int Srch_NZ_Element_Column(RowNum,ColNum)
   { 
   int Count_i;
    {
     for(Count_i=0;Count_i<=GBL_Matrix_Index[RowNum]-1;Count_i++)
       {
        if(ColNum==(GBL_Matrix_ColNos[RowNum][Count_i]))
             {return Count_i;
             }
       }
     return -1;
    }
   }

/***************************************************************************
*
*   FUNCTION    : Locate_NZ_Element_Column()
*   DESCRIPTION : Searches for column number of non-zero element
*                 within each row of A matrix (?).
*
****************************************************************************/

int Locate_NZ_Element_Column(RowNum,ColNum)
     {
      int Count_i;
       {
        for(Count_i=0;Count_i<=FRCD_Matrix_MaxColNos-1;Count_i++)
          {
           if(ColNum==(FORC_Matrix_ColNos[RowNum][Count_i]))
             {return Count_i;
             }
          } 
        return -1;
       }
     }

/***************************************************************************
*
*   FUNCTION    : Determinant_TetraMatrix()
*   DESCRIPTION : Calculates determiment of a 4 x 4 matrix. 
*
****************************************************************************/

double Determinant_TetraMatrix()
{
int a,b,i,j,k,jj,kk,ll,m=0;
double Buff[4][4],Buff1[4],Buff2[4];
double TempDet=0.0,ppp=0.0,pp;
{

for(a=0;a<=3;a++)
 for(b=0;b<=3;b++)
  {Buff[a][b]=Volume_Cord[a][b];}

for(i=0;i<=3;i++) {ll=0;       
 for(j=0;j<=3;j++)  {
  if(j!=i) {

   for(k=0;k<=3;k++)  
    { if((k!=i)&&(k!=j)) {Buff1[m]=k; m++;} }

m=0; jj=Buff1[0];kk=Buff1[1];
pp=Buff[1][j]*(Buff[2][jj]*Buff[3][kk]-Buff[3][jj]*Buff[2][kk]);
ppp+=pp*pow(-1.0,(double)ll); pp=0.0;ll++;
           }
                    }
Buff2[i]=Buff[0][i]*ppp; ppp=0.0;
TempDet+=(Buff2[i])*pow(-1.0,(double)i);
                 }
TempDet=TempDet/6.0; return (TempDet*-1.0);
}
}

/***************************************************************************
*
*   FUNCTION    : swap(x,y)
*   DESCRIPTION : swap the values of x and y
*
****************************************************************************/

void swap(x,y)
        int *x,*y;
        {
        int temp;
        temp = *x; *x=*y; *y=temp;
        }

/***************************************************************************
*   FUNCTION        : EfieldatNode()
****************************************************************************/

void EfieldatNode(X,Y,Z,Out_F)

int X,Y,Z;
FILE *Out_F;
{
  int Count_i,Xcount,Ycount;
  int Xnode, Ynode, Znode;
  Xcount=(int)Dimension_X;Ycount=(int)Dimension_Y;
  Xnode=(X-Min_X);Ynode=(Y-Min_Y);Znode=(Z-Min_Z);
  Count_i=Xnode + ((Xcount+1)*Ynode) + ((Xcount+1)*(Ycount+1)*Znode);
  fprintf(Out_F," %6.3f\t%6.3f\t%6.3f\t",MeshNodeCord[Count_i][0],
  MeshNodeCord[Count_i][1],MeshNodeCord[Count_i][2]);
  fprintf(Out_F," %f\t%f\t%f\t",
  Dat_Buffer[Count_i][0],Dat_Buffer[Count_i][1],Dat_Buffer[Count_i][2]);
  fprintf(Out_F,"\n");

}
/*************************** END OF PROGRAM ********************************/