#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "portmesh.h"
#include "material.h"
#include "bc.h"
#include "tensor.h"
#include "cvtr.h"
#include "Constants.h"
#include "tri.elem"
#include "vtkwriter.h"

portMesh::portMesh(){
  faceCNT = 0;
  edgeCNT = 0;
  nodeCNT = 0;
  objCNT = 0;
  objMAP = 0;
  ndArray = 0;
  fcArray = 0;
  modeIMP = 0.;
  magE = 0.;
  impZ = 0.;
  globToLocMap_ = 0;
}

portMesh::~portMesh(){
  delete [] fcArray;
  delete [] ndArray;
  delete [] objMAP;
  delete globToLocMap_;
}

void portMesh::writeMesh(Material *propMAT, fp_t freq)
{
  FILE *foo;
  char meshName[180];
  int i, j, k, objNum;
  tensor epsr, mur;
  node *nd;
  vtr vv, vv2D;

  sprintf(meshName, "%s.mesh", portName);
  foo = fopen(meshName, "w");
  fprintf(foo, "%.20e %.20e %.20e\n", (coord.O.getx()), (coord.O.gety()), (coord.O.getz()));
  fprintf(foo, "%.20e %.20e %.20e\n", (coord.x_axis.getx()), (coord.x_axis.gety()), (coord.x_axis.getz()));
  fprintf(foo, "%.20e %.20e %.20e\n", (coord.y_axis.getx()), (coord.y_axis.gety()), (coord.y_axis.getz()));
  fprintf(foo, "1.0 \n");
  fprintf(foo, "%d %d %d \n", objCNT, nodeCNT, faceCNT);
  
  for(i = 0; i < nodeCNT; i++)
  {
    nd = ndArray[i];
    vv = nd->getCoord();
    vv2D = coord.Transform(vv);

    fprintf(foo, "0 %.20e %.20e \n", (vv2D.getx()), (vv2D.gety()));
  }

  for(i = 0; i < faceCNT; i++){
    objNum = fcArray[i]->gettetraPtr(0)->getobjNum();
    fprintf(foo, "%d %d %d %d \n", getobjMAP(objNum) + 1, fcArray[i]->getNodePtr(0)->num2D, fcArray[i]->getNodePtr(1)->num2D, fcArray[i]->getNodePtr(2)->num2D);
    for(j = 0; j < 3; j++){
      if(fcArray[i]->getEdge(j) == 0) 
        fprintf(foo, " 1 ");
      else 
        fprintf(foo, " 0 ");
    }

    fprintf(foo, "\n");
  }

  for(i = 0; i < objCNT; i++){
    objNum = objMAP[i];
    epsr = propMAT[objNum].epsr;
    mur = propMAT[objNum].mur;
    epsr = coord.Transform(epsr);

    fprintf(foo, "%d \n", i + 1);
    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++){
        fprintf(foo, "%.20e ", (epsr.getEntry(j, k)));
      }
      fprintf(foo, "\n");
    }

    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++){
        fprintf(foo, "%.20e ", (mur.getEntry(j, k)));
      }
      fprintf(foo, "\n");
    }

    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++){
        fprintf(foo, "0.0 "); 
      }
      fprintf(foo, "\n");
    }

    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++){
        fprintf(foo, "0.0 ");
      }
      fprintf(foo, "\n");
    }
  }

  fprintf(foo, "%d \n", vline.np);
  for(i = 0; i < vline.np; i++){
    vv = vline.coord[i];
    cout << vv.getx() << " " << vv.gety() << " " << vv.getz() << endl;
    vv2D = coord.Transform(vv);
    fprintf(foo, "%.20e %.20e \n", vv2D.getx(), vv2D.gety());
  }
  fclose(foo);
}


void portMesh::writeMesh(Material* propMat){
  int i, j, k;
  vtr vv, vv2D;

  char meshName[STRLEN];
  memset(meshName, 0, STRLEN*sizeof(char));
  sprintf(meshName, "%s.mesh", portName);
  FILE* foo = fopen(meshName, "w");
  vtr& O = coord.getO();
  vtr& xAxis = coord.getXAxis();
  vtr& yAxis = coord.getYAxis();

  fprintf(foo, "%.20e %.20e %.20e\n", (O.getx()), (O.gety()), (O.getz()));
  fprintf(foo, "%.20e %.20e %.20e\n", (xAxis.getx()), (xAxis.gety()), (xAxis.getz()));
  fprintf(foo, "%.20e %.20e %.20e\n", (yAxis.getx()), (yAxis.gety()), (yAxis.getz()));
  fprintf(foo, "1.0\n");
  fprintf(foo, "%d %d %d \n", objCNT, nodeCNT, faceCNT);

  for(i = 0; i < nodeCNT; i++){
    node* nd = ndArray[i];
    vv = nd->getCoord();
    vv2D = coord.Transform(vv);
    fprintf(foo, "0 %.20e %.20e\n", (vv2D.getx()), (vv2D.gety()));
  }

  int objNum = 0;
  for(i = 0; i < faceCNT; i++){
    objNum = fcArray[i]->getTetra(0).getobjNum();
    fprintf(foo, "%d %d %d %d\n", getobjMAP(objNum) + 1,
            globToLocMap_->find(fcArray[i]->getNode(0)->getid())->second,
            globToLocMap_->find(fcArray[i]->getNode(1)->getid())->second,
            globToLocMap_->find(fcArray[i]->getNode(2)->getid())->second);
    for(j = 0; j < 3; j++){
      if(fcArray[i]->getEdge(j)->getbType() == pecType) 
        fprintf(foo, " 1 ");
      else 
        fprintf(foo, " 0 ");
    }
    fprintf(foo, "\n");
  }

  tensor epsr;
  tensor mur;
  for(i = 0; i < objCNT; i++){
    objNum = objMAP[i];
    epsr = propMat[objNum].epsr;
    mur = propMat[objNum].mur;
    epsr = coord.Transform(epsr);

    fprintf(foo, "%d\n", i + 1);
    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++)
        fprintf(foo, "%.20e ", (epsr.getEntry(j, k)));
      fprintf(foo, "\n");
    }
    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++)
        fprintf(foo, "%.20e ", (mur.getEntry(j, k)));
      fprintf(foo, "\n");
    }
    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++)
        fprintf(foo, "%.20e ", 0.0);
      fprintf(foo, "\n");
    }
    for(j = 0; j < 3; j++){
      for(k = 0; k < 3; k++) fprintf(foo, "0.0 ");
      fprintf(foo, "\n");
    }
  }

  fprintf(foo, "%d\n", vline.np);

  for(i = 0; i < vline.np; i++){
    cout << vv.getx() << " " << vv.gety() << " " << vv.getz() << endl;
    vv = vline.coord[i];
    vv2D = coord.Transform(vv);
    fprintf(foo, "%.20e %.20e\n", (vv2D.getx()), (vv2D.gety()));
  }

  fclose(foo);
}

void portMesh::makeObjMap(){
  int i, tmpCNT, objNum;

  tmpCNT = 0;
  for(i = 0; i < faceCNT; i++){
    objNum = fcArray[i]->gettetraPtr(0)->getobjNum();
    tmpCNT = (tmpCNT > objNum) ? tmpCNT : objNum;
  }
  tmpCNT++;

  objMAP = new int[tmpCNT];
  for(i = 0; i < faceCNT; i++){
    objNum = fcArray[i]->gettetraPtr(0)->getobjNum();
    insertOBJ(objNum);
  }
}

int portMesh::getobjMAP(int obj){
  int i;
  for(i = 0; i < objCNT; i++){
    if(objMAP[i] == obj) return i;
  }
  return -1;
}

void portMesh::insertOBJ(int obj){
  int i;

  for(i = 0; i < objCNT; i++){
    if(objMAP[i] == obj) 
      return;
  }

  objMAP[objCNT] = obj;
  objCNT++;
}

void portMesh::setFaceCnt(int faceCnt){
  faceCNT = faceCnt;
  fcArray = new face*[faceCNT];
  memset(fcArray, 0, faceCNT * sizeof(face*));
}

void portMesh::setNodeCnt(int nodeCnt){
  nodeCNT = nodeCnt;
  ndArray = new node*[nodeCNT];
  memset(ndArray, 0, nodeCNT * sizeof(node*));
}

// Read from portnems.vpath the definition of the path which defines
// the voltage integration forthe mode. This is necessary to do 
// de-embedding to 50 ohm transmission lines. ifthe file does not exist
// then the S parameters will be expressed in terms of the intrinsic
// impedances.
//

void portMesh::readVline(fp_t unit_conv){
  FILE *foo;
  char vlineName[180];
  int i, np;
  fp_t x, y, z;
  double x_,y_,z_;

  sprintf(vlineName, "%s.vpath", portName);
  if((foo = fopen(vlineName, "r")) == NULL){
    vline.np = 0;
    vline.coord = 0;

    return;
  }

  fscanf(foo, "%d", &np);
  vline.np = np;
  vline.coord = new vtr[np];
  for(i = 0; i < np; i++){
    fscanf(foo, "%le", &x_);
    fscanf(foo, "%le", &y_);
    fscanf(foo, "%le", &z_);
    (vline.coord[i]).setvtr((fp_t)x_ * unit_conv, (fp_t)y_ * unit_conv, (fp_t)z_ * unit_conv);
  }
  fclose(foo);
}

void portMesh::readVline(fp_t unit_conv, fp_t& xoff, fp_t& yoff, fp_t& zoff){
  FILE *foo;
  char vlineName[180];
  int i, np;
  fp_t x, y, z;
  double x_,y_,z_;

  sprintf(vlineName, "%s.vpath", portName);
  if((foo = fopen(vlineName, "r")) == NULL){
    vline.np = 0;
    vline.coord = 0;

    return;
  }

  fscanf(foo, "%d", &np);
  vline.np = np;
  vline.coord = new vtr[np];
  for(i = 0; i < np; i++){
    fscanf(foo, "%le", &x_);
    fscanf(foo, "%le", &y_);
    fscanf(foo, "%le", &z_);
    (vline.coord[i]).setvtr((fp_t)x_ * unit_conv + xoff,  (fp_t)y_ * unit_conv + yoff,  (fp_t)z_ * unit_conv + zoff);
  }
  fclose(foo);
}

// portMesh::makeRHS
// foreach port, we assemble a right-hand-side to be used as the
// excitation forthe computations of the entire scattering matrix
// of the multiport microwave components.
// void portMesh::makeRHS(fp_t freq, int dim, int *Part, portCounter pC){
//   FILE *foo;
//   char solname[180], cc;
//   int i, j, k, tmpCNT, triMAP[8];
//   vtr  tmpVTR, nvtr[2];
//   vtr hvtr[6], jvtr[6];
//   cVtr tvtr;
//   fp_t jval[8], rhs, cval;
//   face *fc;
//   denseMat elemM(8, 8);
//   Material prop;
//   fp_t scaleF;

//   cval = -Uo;
//   portY.N = dim;
//   portY.entry = new fp_t[dim];

//   sprintf(solname, "%s.1H1", portName);
//   if((foo = fopen(solname, "rb")) == NULL){
//     printf("Could not open File %s \n", solname);
//     exit(1);
//   }
//   fread(&tmpCNT, sizeof(int ), 1, foo);
//   fread(&cc, sizeof(char ), 1, foo);
//   fread(&cc, sizeof(char ), 1, foo);
//   fread(&tmpCNT, sizeof(int ), 1, foo);
//   fread(&gamma, sizeof(fp_t ), 1, foo);
//   fread(&gamma, sizeof(fp_t ), 1, foo);
//   for(i = 0; i < 3; i++){
//     fread(&tmpVTR, sizeof(vtr ), 1, foo);
//   }
//   fread(&modeIMP, sizeof(fp_t ), 1, foo);
//   //  scaleF = modeIMP;
//   //  scaleF = scaleF / impZ;
//   //  scaleF = sqrt(scaleF);
//   scaleF = 1.0;
//   for(i = 0; i < 8; i++) 
//     jval[i] = 0.0;
//   for(i = 0; i < faceCNT; i++){
//     fc = fcArray[i];
//     maketriMAP(fc, Part, triMAP, pC);

//     for(j = 0; j < 6; j++){
//       fread(&tvtr, sizeof(cVtr ), 1, foo);
//       hvtr[j].setvtr(tvtr.x.real(), tvtr.y.real(), tvtr.z.real());
//       hvtr[j] = hvtr[j] * scaleF;

//       jvtr[j].setvtr(hvtr[j].gety() * (-1.0), hvtr[j].getx(), 0.0);
//     }
//     fc->P2toH1(jvtr, jval, coord);
//     fc->makeElemMatrix(prop, coord, 8, &elemM);
//     for(j = 0; j < 8; j++){
//       if(triMAP[j] < 0) 
//         continue;
//       rhs = 0.0;
//       for(k = 0; k < 8; k++){
//         rhs = rhs + jval[k] * elemM.getEntry(j, k);
//       } 
//       portY.addentry(triMAP[j], rhs * cval);
//     }
//   }
//   fclose(foo);
// }

struct dbl_vtr{
  double x,y,z;
};

struct dbl_cvtr{
  complex<double> x,y,z;
};



// -----------------------------------------------------------------------------
// portMesh::makeRHS_E
// Reads electric-field port samples from "<portName>.1E1", one block per face,
// 6 complex vector samples per face (second-order triangle).
//
// Each sample dbl_cvtr has .x/.y/.z as std::complex<double>.
// Samples are rotated into global coordinates via coord.x_axis/y_axis/z_axis
// and scaled by scaleF (currently 1).
// -----------------------------------------------------------------------------
void portMesh::makeRHS_E()
{
  FILE *foo;
  char solname[StrOutput], cc;
  int i, j, tmpCNT, triMAP[8];
  vtr nvtr[2];
  dbl_vtr tmpVTR;
  vtr evtr[NumOfEdges];  // holds 6 E samples for the current face
  dbl_cvtr tvtr;
  fp_t cval;
  face *fc;
  fp_t scaleF;

  double gamma_,modeIMP_;

  sprintf(solname, "%s.1E1", portName);

  cout << "Reading port mode solution from file: " << solname << endl;
  
  
  if((foo = fopen(solname, "rb")) == NULL)
  {
    printf("Could not open File %s \n", solname);
    exit(1);
  }


  fread(&tmpCNT, sizeof(int), 1, foo);
  fread(&cc, sizeof(char), 1, foo);
  fread(&cc, sizeof(char), 1, foo);
  fread(&tmpCNT, sizeof(int), 1, foo);
  fread(&gamma_, sizeof(double), 1, foo);
  fread(&gamma_, sizeof(double), 1, foo);
  for(i = 0; i < NumOfUnitaryVectors; i++)
  {
    fread(&tmpVTR, sizeof(dbl_vtr), 1, foo);
  }
  fread(&modeIMP_, sizeof(double), 1, foo);

  gamma = (fp_t)gamma_;
  modeIMP = (fp_t)modeIMP_;

  // --- per-face blocks ---
  scaleF = 1.0;
  for(j = 0; j < NumOfEdges; j++)
  {
    evtr[j].setvtr(0.0, 0.0, 0.0);
  }

  // read 6 complex-vector samples (x,y,z) for this face
  for(i = 0; i < faceCNT; i++)
  {
    fc = fcArray[i];
    for(j = 0; j < 6; j++)
    {
      fread(&tvtr, sizeof(dbl_cvtr), 1, foo);
      
      // rotate from file-local axes to global axes
      evtr[j] = coord.x_axis * (fp_t)(tvtr.x.real()) + coord.y_axis * (fp_t)(tvtr.y.real()) + coord.z_axis * (fp_t)(tvtr.z.real());
      evtr[j] = evtr[j] * scaleF;
    }

    // convert the 6 P2 samples to H1 coefficients and store on the face
    fc->Set_WgSolverExcitation_E(evtr);
  }
  fclose(foo);
}



void portMesh::makeRHS_H()
{
  FILE *foo;
  char solname[180], cc;
  int i, j, tmpCNT, triMAP[8];
  vtr nvtr[2];
  dbl_vtr tmpVTR;
  vtr hvtr[6];
  dbl_cvtr tvtr;
  fp_t cval;
  face *fc;
  fp_t scaleF;
  double gamma_,modeIMP_;

  sprintf(solname, "%s.1H1", portName);
  if((foo = fopen(solname, "rb")) == NULL){
    printf("Could not open File %s \n", solname);
    exit(1);
  }
  fread(&tmpCNT, sizeof(int), 1, foo);
  fread(&cc, sizeof(char), 1, foo);
  fread(&cc, sizeof(char), 1, foo);
  fread(&tmpCNT, sizeof(int), 1, foo);
  fread(&gamma_, sizeof(double), 1, foo);
  fread(&gamma_, sizeof(double), 1, foo);

  for(i = 0; i < 3; i++)
  {
    fread(&tmpVTR, sizeof(dbl_vtr), 1, foo);
  }
  
  fread(&modeIMP_, sizeof(double ), 1, foo);

  gamma = (fp_t)gamma_;
  modeIMP = (fp_t)modeIMP_;

  scaleF = 1.0;
  for(j = 0; j < 6; j++){
    hvtr[j].setvtr(0.0, 0.0, 0.0);
  }
  for(i = 0; i < faceCNT; i++){
    fc = fcArray[i];
    for(j = 0; j < 6; j++){
      fread(&tvtr, sizeof(dbl_cvtr), 1, foo);
      hvtr[j] = coord.x_axis * (fp_t)(tvtr.x.real()) + coord.y_axis * (fp_t)(tvtr.y.real()) + coord.z_axis * (fp_t)(tvtr.z.real());
      hvtr[j] = hvtr[j] * scaleF;
    }
    fc->Set_WgSolverExcitation_H(hvtr);
  }
  fclose(foo);
}

void portMesh::makeCoordSystem()
{
  vtr zAxis = findFaceNormal(*(fcArray[0]));
  zAxis = zAxis * (-1.0); //point inside tetr(0) 
  // zAxis = zAxis *(1.0); // forSMT
  // forthe SMT example somehow the zAxis was (0,0,1) but from the geometry 
  // it should be (0,0,-1). So I flip the sign to make the right direction 
  coord.setz_axis(zAxis);   // I always assume the zAxis is the portDirection

  node* vertexNode = 0;
  node* edgeNode = 0;
  findNodes(&vertexNode, &edgeNode);
  if(vertexNode == 0){
    cerr << "ERROR: PortMesh::makeCoordSystem: " << "Could not find a vertex node forport " << portName << endl;
    exit(1);
  }
  if(edgeNode == 0){
    cerr << "ERROR: PortMesh::makeCoordSystem: " << "Could not find an edge node forport " << portName << endl;
    exit(1);
  }
  vtr O = vertexNode->getCoord();
  coord.setO(O);

  vtr xAxis = edgeNode->getCoord() - O;
  coord.setx_axis(xAxis);

  vtr yAxis = zAxis * xAxis;
  coord.sety_axis(yAxis);
  vtr PortDirection = coord.getZAxis();

  cout << portName << endl;
  cout << "( " << PortDirection.getx() << ", " << PortDirection.gety() << ", " << PortDirection.getz() << " )" << endl;
  
  fp_t Port_Direction_mag = sqrt(PortDirection.getx() * PortDirection.getx() + PortDirection.gety() * PortDirection.gety() + PortDirection.getz() * PortDirection.getz());

  fp_t px = PortDirection.getx()/Port_Direction_mag;
  fp_t py = PortDirection.gety()/Port_Direction_mag;
  fp_t pz = PortDirection.getz()/Port_Direction_mag;
  PortDirection_vtr.setvtr(px,py,pz);

  for(int i = 0; i < faceCNT; i++) 
    fcArray[i]->getbcPtr()->set_PortDirection( PortDirection_vtr.getx(), PortDirection_vtr.gety(), PortDirection_vtr.getz());
}

vtr portMesh::findFaceNormal(face& Face){
  int i;

  vtr vv0 = Face.getNode(0)->getCoord();
  vtr vv1 = Face.getNode(1)->getCoord();
  vtr vv2 = Face.getNode(2)->getCoord();
  vtr nvtr = (vv1 - vv0) * (vv2 - vv0);
  nvtr = nvtr * 0.5;

  tetra& Tetra = Face.getTetra(0);
  if(Tetra == 0){
    cerr << "ERROR: DgGroup::findFaceNormal: " << "Port face does not have a correct tetra link." << endl;
    exit(1);
  }else{
    vv0.setvtr(0., 0., 0.);
    for(i = 0; i < 4; i++){
      vv0 = vv0 + Tetra.getNode(i)->getCoord();
    }
    vv0 = vv0 * 0.25;
    vv0 = vv0 - vv1;
    if(dotP(vv0, nvtr) > 0.) 
      nvtr = nvtr * (-1.);
  }
  nvtr = nvtr * (1.0 / nvtr.magnitude());
  return nvtr;
}

void portMesh::findNodes(node** vNode, node** eNode){
  int i, j;
  vtr vv1;
  vtr vv2;

  for(i = 0; i < faceCNT; i++){
    bool foundVertex = false;
    bool foundEdge = false;
    node* vertexNode = 0;
    node* edgeNode = 0;
    for(j = 0; j < 3; j++){
      node* Node = fcArray[i]->getNodePtr(j);
      if(Node->getPType() == VERTEX_NODE){
        if(foundVertex == true){
          if(foundEdge == true){
            vv1 = edgeNode->getCoord() - vertexNode->getCoord();
            vv2 = Node->getCoord() - vertexNode->getCoord();
            if(vv2.magnitude() < vv1.magnitude()){
              edgeNode = Node;
            }
          }else{
            foundEdge = true;
            edgeNode = Node;
          }
        }else{
          vertexNode = Node;
          foundVertex = true;
        }
      }
      if(Node->getPType() == EDGE_NODE){
        foundEdge = 1;
        edgeNode = Node;
      }
    }
    if((foundVertex) && (foundEdge)){
      *vNode = vertexNode;
      *eNode = edgeNode;
      return;
    }
  }
}

void portMesh::writeVtk()
{
  int i, j;
  char  hFileName[STRLEN];
  char  eFileName[STRLEN];
  FILE* fdE;
  FILE* fdH;
  char cc;
  int dummy;
  dbl_vtr dummyVtr;
  complex<double> gamma_;
  double modeIMP_;

  sprintf(eFileName, "%s.1E1", portName);
  if((fdE = fopen(eFileName, "rb")) == 0)
  {
    cerr << "ERROR: PortMesh::makeRhs: " << "Could not open file " << eFileName << " forread." << endl;
    exit(1);
  }
  
  fread(&dummy, sizeof(int), 1, fdE);
  fread(&cc, sizeof(char), 1, fdE);
  fread(&cc, sizeof(char), 1, fdE);
  fread(&dummy, sizeof(int), 1, fdE);
  fread(&gamma_, sizeof(complex<double>), 1, fdE);
  for(i = 0; i < 3; i++) 
    fread(&dummyVtr, sizeof(dbl_vtr), 1, fdE);
  fread(&modeIMP_, sizeof(double), 1, fdE);

  sprintf(hFileName, "%s.1H1", portName);
  if((fdH = fopen(hFileName, "rb")) == 0)
  {
    cerr << "ERROR: PortMesh::makeRhs: " << "Could not open file " << hFileName << " forread." << endl;
    exit(1);
  }
  
  fread(&dummy, sizeof(int), 1, fdH);
  fread(&cc, sizeof(char), 1, fdH);
  fread(&cc, sizeof(char), 1, fdH);
  fread(&dummy, sizeof(int), 1, fdH);
  fread(&gamma_, sizeof(complex<double>), 1, fdH);
  for(i = 0; i < 3; i++) 
    fread(&dummyVtr, sizeof(dbl_vtr), 1, fdH);
  fread(&modeIMP_, sizeof(double), 1, fdH);

  // fill the port field with averaged values
  vtr* eField = new vtr[nodeCNT];
  vtr* hField = new vtr[nodeCNT];
  int* count = new int[nodeCNT];
  memset(count, 0, nodeCNT*sizeof(int));

  int index;
  dbl_cvtr evtr;
  dbl_cvtr hvtr;
  vtr  tmpE; // change here
  vtr  tmpH;
  for(i = 0; i < faceCNT; i++){
    for(j = 0; j < SEC_ORDER_FACE_NODES; j++){
      fread(&evtr, sizeof(dbl_cvtr), 1, fdE);
      fread(&hvtr, sizeof(dbl_cvtr), 1, fdH);

      if(j < 3)
      {
        index = globToLocMap_->find(fcArray[i]->getNode(j)->getid())->second;
        
        tmpE = coord.x_axis * (fp_t)(evtr.x.real()) + coord.y_axis * (fp_t)(evtr.y.real()) + coord.z_axis * (fp_t)(evtr.z.real());
        tmpH = coord.x_axis * (fp_t)(hvtr.x.real()) + coord.y_axis * (fp_t)(hvtr.y.real()) + coord.z_axis * (fp_t)(hvtr.z.real());
        
        eField[index] = eField[index] + tmpE;
        hField[index] = hField[index] + tmpH;
        (count[index])++;
      }
    }
  }

  fclose(fdE);
  fclose(fdH);

  for(i = 0; i < nodeCNT; i++){
    eField[i] = eField[i] / static_cast<fp_t>(count[i]);
    hField[i] = hField[i] / static_cast<fp_t>(count[i]);
  }

  node** locNodeArray = new node*[nodeCNT];
  for(i = 0; i < nodeCNT; i++){
    node& Node = *(ndArray[i]);
    int index = globToLocMap_->find(Node.getid())->second;
    locNodeArray[index] = new node(index, Node.getPType(), Node.getSingOrder(), Node.getCoord().getx(), Node.getCoord().gety(), Node.getCoord().getz());
  }

  face** locFaceArray = new face*[faceCNT];
  for(i = 0; i < faceCNT; i++){
    face& Face = *(fcArray[i]);
    locFaceArray[i] = new face(Face);
    locFaceArray[i]->setFace(
        locNodeArray[globToLocMap_->find(Face.getNode(0)->getid())->second],
        locNodeArray[globToLocMap_->find(Face.getNode(1)->getid())->second],
        locNodeArray[globToLocMap_->find(Face.getNode(2)->getid())->second]);
  }

  VtkWriter vtkWriter(1.);
  vtkWriter.writeTriUg(portName, nodeCNT, locNodeArray, faceCNT, locFaceArray, eField, hField, 1);

  for(i = 0; i < nodeCNT; i++) 
    delete locNodeArray[i];
  delete [] locNodeArray;
  for(i = 0; i < faceCNT; i++) 
    delete locFaceArray[i];
  delete [] locFaceArray;
  delete [] eField;
  delete [] hField;
  delete [] count;
}




// -----------------------------------------------------------------------------
// portMesh::makeRHS_TEM  (HFSS-style symmetry port)
//
// Builds a uniform TEM-like snapshot on the port cross-section:
//
//   • Symmetry planes are applied only on the port boundary (not physical walls).
//   • E is set to point across the PEC symmetry plane (normal to it, in-plane).
//   • H = (1/η) (ez × E),   with ez the port normal (power +S).
//   • |E| is fixed to 1 V/m  (no 1-W renormalization).
//
// Inputs:
//   freq_Hz, eps_r : frequency [Hz] and relative permittivity (uniform medium)
//   Ex_dir_*       : a vector pointing normal to the PEC symmetry plane,
//                    lying in the port plane (if not, it will be projected).
//   PortDir_*      : the port normal (+S direction).
//
// Side effects / state:
//   • For every face, provides 6 samples (2nd-order tri) via Set_WgSolverExcitation_E/H.
//   • Stores wave (field) impedance η = sqrt(μ/ε) in modeIMP and modeImpedance.
//   • Stores gamma = k = ω sqrt(μ ε) (lossless).
//   • Logs the resulting power for |E|=1 V/m, i.e.  P = 0.5 * (1/η) * A.
// -----------------------------------------------------------------------------
void portMesh::makeRHS_TEM(double freq_Hz, double eps_r,
                           double Ex_dir_x, double Ex_dir_y, double Ex_dir_z,   // E direction (normal to PEC plane, in-plane)
                           double PortDir_x, double PortDir_y, double PortDir_z)// propagation (+S)
{
  using std::complex;

  std::cout << "  [TEM|Sym] Fill for " << portName
            << " @ f=" << freq_Hz/1e9 << " GHz, eps_r=" << eps_r
            << "  (|E| = 1 V/m)\n";

  // --- 1) Physical constants and medium parameters (SI) ----------------------
  const double c0   = 299792458.0;
  const double mu0  = 4e-7 * M_PI;
  const double eps0 = 1.0 / (mu0 * c0 * c0);

  const double mu   = mu0;                 // non-magnetic medium
  const double eps  = eps_r * eps0;
  const double w    = 2.0 * M_PI * freq_Hz;
  const double k    = w * std::sqrt(mu * eps);   // β magnitude for lossless TEM
  const double eta  = std::sqrt(mu / eps);       // wave impedance E/H (field impedance)

  // Store in your legacy members
  gamma          = static_cast<fp_t>(k);
  modeIMP        = static_cast<fp_t>(eta);
  modeImpedance  = complex<double>(eta, 0.0);

  // --- 2) Local orthonormal frame {ex, ey, ez} --------------------------------
  // ez = port normal (+S direction)
  vtr ez(PortDir_x, PortDir_y, PortDir_z);
  ez = normalize(ez);
  if (norm(ez) == 0.0) {
    std::cerr << "  [TEM|Sym] ERROR: PortDir is zero.\n";
    return;
  }

  // ex = desired E direction (normal to PEC symmetry plane), projected to port plane
  vtr ex(Ex_dir_x, Ex_dir_y, Ex_dir_z);
  if (norm(ex) == 0.0) {
    std::cerr << "  [TEM|Sym] ERROR: E direction is zero.\n";
    return;
  }
  // Project onto plane ⟂ ez and normalize
  ex = ex - ez * dot(ex, ez);
  ex = normalize(ex);
  if (norm(ex) == 0.0) 
  {
    // Given direction was (nearly) parallel to ez; pick a stable in-plane axis
    std::cerr << "  [TEM|Sym] WARNING: E direction ~ parallel to PortDir; choosing in-plane axis.\n";
    vtr up = (std::fabs(ez.getz()) > 0.9) ? vtr(0,1,0) : vtr(0,0,1);
    ex = normalize( cross(up, ez) );
    if (norm(ex) == 0.0) { up = vtr(0,1,0); ex = normalize( cross(up, ez) ); }
    if (norm(ex) == 0.0) {
      std::cerr << "  [TEM|Sym] ERROR: cannot establish in-plane E direction.\n";
      return;
    }
  }

  // ey completes right-handed frame
  vtr ey = cross(ez, ex);

  // --- 3) Compute port area A (sum of triangle areas) -------------------------
  // OPTION A: use node* (NOT const node*) to call non-const node::getCoord()
  double A_port = 0.0;
  for (int f = 0; f < faceCNT; ++f) {
    node* n0 = fcArray[f]->getNode(0);
    node* n1 = fcArray[f]->getNode(1);
    node* n2 = fcArray[f]->getNode(2);

    const vtr r0(n0->getCoord().getx(), n0->getCoord().gety(), n0->getCoord().getz());
    const vtr r1(n1->getCoord().getx(), n1->getCoord().gety(), n1->getCoord().getz());
    const vtr r2(n2->getCoord().getx(), n2->getCoord().gety(), n2->getCoord().getz());

    A_port += tri_area(r0, r1, r2);
  }
  if (!(A_port > 0.0)) {
    std::cerr << "  [TEM|Sym] ERROR: port area is zero.\n";
    return;
  }

  // --- 4) Fix |E| = 1 V/m; H consistent with η; compute resulting power -------
  const double E0  = 1.0;                        // requested fixed magnitude
  const vtr    Evec = ex * static_cast<fp_t>(E0);
  const vtr    Hvec = ey * static_cast<fp_t>(E0 / eta);

  // For info: total time-average power P = 0.5 * (E×H)·ez * A
  // Here (E×H)·ez = (E0^2)/η and is uniform over the port
  const double P_result = 0.5 * ( (E0 * E0) / eta ) * A_port;

  // --- 5) Provide 6 samples per face (uniform, so same vector at each sample) -
  vtr E6[6], H6[6];
  for (int f = 0; f < faceCNT; ++f) 
  {
    E6[0] = Evec;  H6[0] = Hvec;  // vertex 0
    E6[1] = Evec;  H6[1] = Hvec;  // vertex 1
    E6[2] = Evec;  H6[2] = Hvec;  // vertex 2
    E6[3] = Evec;  H6[3] = Hvec;  // midside (0,1)
    E6[4] = Evec;  H6[4] = Hvec;  // midside (1,2)
    E6[5] = Evec;  H6[5] = Hvec;  // midside (2,0)

    fcArray[f]->Set_WgSolverExcitation_E(E6);
    fcArray[f]->Set_WgSolverExcitation_H(H6);

    cout << "    Face " << f << ": E = ("
         << Evec.getx() << ", " << Evec.gety() << ", " << Evec.getz() << ") V/m,  |E|=" << norm(Evec) << " V/m\n";
  }

  // --- 6) Log summary ---------------------------------------------------------
  std::cout << "  [TEM|Sym] η = " << eta
            << " Ω (E/H),  A = " << A_port << " m^2,  |E| = " << E0 << " V/m"
            << ",  resulting P = " << P_result << " W\n";
}