#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <cmath>
#include <cstring>
#include <ctime>
#include <iostream>
#include <fstream>

#include <unistd.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <sys/times.h>

#include "global.h"
#include "proc.h"
#include "Constants.h"
#include "float_regcvtr.h"
#include "group.h"
#include "MultiLevelSVD.h"
#include "n2f_dd.h"
#include "transform.h"
#include "integral.h"
using namespace std;

// global variables
// command line parameters
#ifndef ZERO_7
double const static ZERO_7 = 1.e-7;
#endif
vtr O;

double unitConv;
tri3d **groupTriArray;
regcvtr *groupRegJ, *groupRegM;
fftwf_complex *workspace, *field, *out;
fftwf_complex *tmp_x, *tmp_y, *tmp_z;
cVtr *Epart;

double MagE( Complex eTheta, Complex ePhi )
{
  double mag = norm(eTheta) + norm(ePhi);

  return sqrt(mag);
}

int main(int argc, char *argv[])
{
  tms tmsStart, tmsEnd;
  times(&tmsStart);

  int i, j;
  int numTheta, numPhi,OPT,MODE;
  char prjName[180];
  double freq;
  
  // Read in input parameters
 if ( argc == 10 ) {
    strcpy( prjName, argv[1] );
    freq = atof( argv[2] );
    O.setvtr( atof(argv[3]), atof(argv[4]), atof(argv[5]) );
    numTheta = atoi( argv[6] );
    numPhi = atoi( argv[7] );
    OPT = atoi(argv[8]);
    MODE=atoi(argv[9]);
  } else {
    cout << "Input Project Name: ";
    cin >> prjName;
    cout << "Input Frequency (in MHz): ";
    cin >> freq;
    double Ox, Oy, Oz;
    cout << "Input Coordinate System Center X (meters): ";
    cin >> Ox;
    cout << "Input Coordinate System Center Y (meters): ";
    cin >> Oy;
    cout << "Input Coordinate System Center Z (meters): ";
    cin >> Oz;
    O.setvtr( Ox, Oy, Oz );
    cout << "Input Number of Sample Points along Theta: ";
    cin >> numTheta;
    cout << "Input Number of Sample Points along Phi: ";
    cin >> numPhi;
    cout << "Option for file format: 0(magnitude), 1(real and imaginary) ";
    cin >> OPT;
    cout<<"Input the mode:0(radiation),1(scattering)";
    cin>>MODE;
  }

  // scattering
  int n;
  double EincMag=1.;
  double theta, phi, thetaInc, phiInc, theta_in_rad, phi_in_rad;
  cVtr Einc, Hinc, einc, hinc;
  vtr k;
  double ko = TwoPi * freq * MEGA / Vo;
  double deltaTheta = 180.0 / double(numTheta);
  double deltaPhi = 360.0 / double(numPhi);
  Complex **eTheta, **ePhi;
  eTheta = new Complex *[numTheta + 1];
  ePhi = new Complex *[numTheta + 1];
  for (i = 0; i < numTheta + 1; i ++) {
    eTheta[i] = new Complex [numPhi + 1];
    ePhi[i] = new Complex [numPhi + 1];
  }


  vtr thetaVtr, phiVtr;

  int regions = 0, *FEM, *ID;
  string* 	regionName;
  ReadRegions(prjName, &FEM, &ID, regions,&regionName);

 cout<<"regions"<<regions<<endl;

  bool Pec[regions];

//   for (i = 0; i < regions; i ++) {
//     Pec[i] = false;
//     if (FEM[i] == 0 || FEM[i] == 2) Pec[i] = true;
//   }

     bool SecondOrder[regions];

        for (i = 0; i < regions; i ++) {
                Pec[i] = false;
                SecondOrder[i]=true;
                cout<<"Region "<<i<<"is: ";
                if (FEM[i] == 0 || FEM[i] == 2){
                        Pec[i] = true;
                        SecondOrder[i]=false;
                        cout<<Pec[i]<<"MOM\n";
                }
                else if (FEM[i] == 3){
                        SecondOrder[i]=false;
                        cout<<"FEBI\n";
                }
                else{
                        cout<<"FEM\n";
                }
        }

   // read tri and currents
  int numAbcTri[regions];
  tri3d***	abcTriArray = new tri3d **[regions];
  regcvtr**	regJ = new regcvtr *[regions];
  regcvtr**	regM = new regcvtr *[regions];
  char grpName[180], meshName[180], JcurName[180], McurName[180];


  int np;
  for (i = 0; i < regions; i ++) {

//      cout<<"PEC[i] 0"<<Pec[i]<<endl;

//     if (ID[i] == 0) sprintf(grpName, "%s", prjName);
//     else sprintf(grpName, "%s%d", prjName, ID[i]);
	  sprintf(grpName, "%s", regionName[i].c_str());
    np = 6;
//    if (Pec[i]) np = 3; // PO first order
//     if (FEM[i] == 0) np = 3; // PO second order

     if (SecondOrder[i]==false) np = 3;

//      cout<<"PEC[i] 1"<<Pec[i]<<endl;

    sprintf(meshName, "%s", grpName);
    makeTriArray(numAbcTri[i], &(abcTriArray[i]), meshName);
    sprintf(JcurName, "%s.curJ", grpName);
    setRegister(&(regJ[i]), numAbcTri[i], abcTriArray[i], JcurName, np);

//      cout<<"PEC[i] 2"<<Pec[i]<<endl;

    if (!Pec[i]) {
      sprintf(McurName, "%s.curM", grpName);
      setRegister(&(regM[i]), numAbcTri[i], abcTriArray[i], McurName, np);
    }else
	    regM[i]=NULL;
  }

  InitFFTW();

  // partition the meshing
  group<tri3d> **top = new group<tri3d> *[regions];
  group<tri3d> ***SVDgroup = new group<tri3d> **[regions];
  int NumSVDgroup[regions], grpSize[regions];
  for (i = 0; i < regions; i ++) {
    top[i] = InitilizeTopGroup( abcTriArray[i], numAbcTri[i] );
    grpSize[i] = int(1.5 * sqrt(double(numAbcTri[i])));
    SVDgroup[i] = MeshPartition(top[i], grpSize[i], NumSVDgroup[i]);
  }

  int maxTriNum = 0;
  int maxSizePhi = 0, maxSizeTheta = 0;
  int sizeTheta, sizePhi;
  for (i = 0; i < regions; i ++) {
    for ( j = 0; j < NumSVDgroup[i]; j ++ ) {
      int nTri = SVDgroup[i][j]->GetnElement();
      if (nTri > maxTriNum) maxTriNum = nTri;
      SVDgroup[i][j]->ComputeGroupRadius();
      double d = 2.0 * SVDgroup[i][j]->GetRadius();
      int Lp = int( ko * d + 5.0 * log(ko * d + Pi) );
      sizeTheta = 2 * Lp + 1;
      sizePhi = 2 * Lp;
      sizePhi = (sizePhi > numPhi) ? numPhi : sizePhi;
      sizeTheta = (Lp > numTheta) ? 2*numTheta  : sizeTheta;
      if (sizePhi > maxSizePhi) maxSizePhi = sizePhi;
      if (sizeTheta > maxSizeTheta) maxSizeTheta = sizeTheta;
    }
  }

  groupTriArray = new tri3d*[maxTriNum];

  groupRegJ = new regcvtr [maxTriNum];
  groupRegM = new regcvtr [maxTriNum];
  for (i = 0; i < maxTriNum; i ++) {
    groupRegJ[i].init(6);
    groupRegM[i].init(6);
  }

  int coeffSize = maxSizeTheta * maxSizePhi;
  Epart = new cVtr [coeffSize];
  tmp_x = new fftwf_complex [coeffSize];
  tmp_y = new fftwf_complex [coeffSize];
  tmp_z = new fftwf_complex [coeffSize];

  // allocate storage space 
  sizeTheta = 2 * numTheta;
  sizePhi = numPhi;
  coeffSize = sizeTheta * sizePhi;
  out = new fftwf_complex [coeffSize];
  field = new fftwf_complex [coeffSize];
  workspace = new fftwf_complex [coeffSize];

  cVtr **MainField = new cVtr *[regions];
  for (i = 0; i < regions; i ++) {

    MainField[i] = new cVtr[coeffSize];

    // calculate far field
    for ( j = 0; j < NumSVDgroup[i]; j ++ ) {
      if (j % 5 == 0) 
	cout << "group = " << j << " / " << NumSVDgroup[i] << endl;
      FastFarField( SVDgroup[i], j, ko, regJ[i], regM[i], numTheta, numPhi, 
		    MainField[i], Pec[i] );
    }
  }
  delete [] out; delete [] field; delete [] workspace; delete [] groupTriArray;
  delete [] Epart; delete [] tmp_x; delete [] tmp_y; delete [] tmp_z;
  for (i = 0; i < maxTriNum; i ++) {
    groupRegJ[i].del();
    groupRegM[i].del();
  }
  delete [] groupRegJ; delete [] groupRegM;

  // split far field into theta and phi components
  int idx = 0;
  cVtr E;
  for ( i = 0; i <= numTheta; i ++ ) {
    theta = double(i) * deltaTheta;
    for ( j = 0; j < numPhi; j ++ ) {
      phi = double(j) * deltaPhi;

      thetaVtr = ThetaDirection(theta, phi);
      phiVtr = PhiDirection(theta, phi);
      E.reset();
      for (n = 0; n < regions; n ++) E = E + MainField[n][idx];
      eTheta[i][j] = dotP( E, thetaVtr );
      ePhi[i][j] = dotP( E, phiVtr );

      idx ++;
    }
    eTheta[i][numPhi] = eTheta[i][0];
    ePhi[i][numPhi] = ePhi[i][0];
  }

    for (i = 0; i < regions; i ++)
          cout<< MainField[i][0].getx() << MainField[i][1].getx()<<endl;

  for (i = 0; i < regions; i ++) delete [] MainField[i];
  delete [] MainField;

if(MODE==0){  

double Gmax = 1.0;

//  for (i = 0; i <= numTheta; i ++) {
//    for (j = 0; j <= numPhi; j ++) {
//      double eMag = MagE(eTheta[i][j], ePhi[i][j]);
	
//      Gmax = (eMag > Gmax) ? eMag : Gmax;
//    }
//  }
//  for (i = 0; i <= numTheta; i ++) {
//    for (j = 0; j <= numPhi; j ++) {
//      eTheta[i][j] /= Gmax;
//      ePhi[i][j] /= Gmax;
//    }
//  }

  // Compute the beam solid angle
  double omegaA = 0.0;
//    for (i = 0; i < numTheta; i ++) {
//      double theta_in_rad = (deltaTheta * (double)(i) + 0.5 * deltaTheta) * 
//        Pi / 180.0;
//     
//      for (j = 0; j < numPhi; j ++) {
//        Complex thetaValue = (eTheta[i][j] + eTheta[i][j + 1] + 
//  			    eTheta[i + 1][j] + eTheta[i + 1][j + 1]) * 0.25f;
//        Complex phiValue = (ePhi[i][j] + ePhi[i][j + 1] + ePhi[i + 1][j] +
//  			  ePhi[i + 1][j + 1]) * 0.25f;
// 
//        omegaA += (norm(thetaValue) + norm(phiValue)) * sin(theta_in_rad);
//      }
//    }

 for (i = 0; i < numTheta; i ++) {
   double theta_in_rad = deltaTheta * (double)(i) * Pi / 180.0;
    
   for (j = 0; j < numPhi; j ++) {
     Complex thetaValue = eTheta[i][j];
     Complex phiValue = ePhi[i][j];
      
     omegaA += (norm(thetaValue) + norm(phiValue)) * sin(theta_in_rad);
   }
 }

  omegaA = omegaA * ( deltaTheta * Pi / 180.0) * (deltaPhi * Pi / 180.0);

  double Prad = omegaA * Gmax * Gmax / No / 2.0;

//  omegaA = 4098.68;
//  double Prad = 5.45678;

  cout << "Omega " << omegaA << endl;
  cout<< "No "<< No <<endl;
   cout << "total power = " << Prad << endl;
    double Pavg = Prad / 4.0 / Pi;
    cout << "average power = " << Pavg << endl;

  double directivity = 4.0 * Pi / omegaA;
  cout << "directivity = " << 10 * log10(directivity) << " dB" << endl;

  char antName[180];
  sprintf(antName, "%s.ap", prjName);
  ofstream fap(antName, ios::out);
  fap << freq << " " << numTheta << " " << numPhi << endl;

  for (i = 0; i <= numTheta; i ++) {
    double theta = deltaTheta * (double)(i);
    for (j = 0; j <= numPhi; j ++) {
      double phi = deltaPhi * (double)(j);

//      eTheta[i][j] *= sqrt(directivity);
//      ePhi[i][j]   *= sqrt(directivity);

      if (OPT == 0) {
	fap << theta << " " << phi << " " 
//	    << norm(eTheta[i][j]) << " " << norm(ePhi[i][j]) << endl;
//            << directivity * norm(eTheta[i][j]) << " "
//            << directivity * norm(ePhi[i][j])   << endl;

            << eTheta[i][j].real() << " "<< (eTheta[i][j]).imag() << " "
            << ePhi[i][j].real()   << " " << ePhi[i][j].imag() << " "
            << (norm(eTheta[i][j])+norm(ePhi[i][j])) / 2.0 / No / Pavg<<endl;


      } else {
	fap << theta << " " << phi << " " 
	    << (eTheta[i][j]).real() << " " << (eTheta[i][j]).imag() << " "
	    << (ePhi[i][j]).real()   << " " << (ePhi[i][j]).imag()   << endl;
      }
    }
  }

  fap.close();

    sprintf(antName, "%s.dBi", prjName );
    fap.open(antName, ios::out);
    fap << freq << " " << numTheta << " " << numPhi << endl;

    for (i = 0; i <= numTheta; i++ ) {
      theta = deltaTheta * double(i);
      for (j = 0; j <= numPhi; j++ ) {
        phi = deltaPhi * double(j);
        fap << theta << " " << phi << " "
            << norm(eTheta[i][j]) * Gmax * Gmax / Pavg << " "
            << norm(ePhi[i][j]) * Gmax * Gmax / Pavg << " "
            << (norm(eTheta[i][j])+norm(ePhi[i][j])) * Gmax * Gmax/Pavg<<endl;
      }
    }
    fap.close();


}else if (MODE==1){
double sigma_th_th, sigma_ph_th, sigma_th_ph, sigma_ph_ph;
    double sigma_th_th_dbsw, sigma_ph_th_dbsw;
    double sigma_th_ph_dbsw, sigma_ph_ph_dbsw;
    double sigma, sigma_dbsw;
    double sigma_th, sigma_ph;

    char scatName[360];
    sprintf(scatName, "%s.cs", prjName);
    FILE *fsc = fopen(scatName, "wt");
    fprintf(fsc, "%5.3f %d %d\n" , freq, numTheta, numPhi );

    for (i = 0; i <= numTheta; i ++) {
      theta = double(i)*deltaTheta;

      for (j = 0; j <= numPhi; j ++) {
	phi = double(j)*deltaPhi;

	double eThMag = sqrt(norm(eTheta[i][j]));
	double ePhMag = sqrt(norm(ePhi[i][j]));
	double eMag = MagE(eTheta[i][j], ePhi[i][j]);

	thetaVtr = ThetaDirection(theta, phi);
	phiVtr = PhiDirection(theta, phi);

	Complex eThetaInc = dotP(Einc, thetaVtr);
	Complex ePhiInc = dotP(Einc, phiVtr);
	double EincThMag = sqrt(norm(eThetaInc));//EincMag
	double EincPhMag = sqrt(norm(ePhiInc));

 	sigma_th_th = (eThMag * eThMag)/(EincMag*EincMag)/(4.*Pi);
 	sigma_ph_th = (ePhMag * ePhMag)/(EincMag*EincMag)/(4.*Pi);
 	sigma_th_ph = (eThMag * eThMag)/(EincMag*EincMag)/(4.*Pi);
 	sigma_ph_ph = (ePhMag * ePhMag)/(EincMag*EincMag)/(4.*Pi);
 
 	sigma_th_th_dbsw = 10.*log10(sigma_th_th);
 	sigma_ph_th_dbsw = 10.*log10(sigma_ph_th);
 	sigma_th_ph_dbsw = 10.*log10(sigma_th_ph);
 	sigma_ph_ph_dbsw = 10.*log10(sigma_ph_ph);
 	sigma = 4.*Pi * (eMag * eMag)/(EincMag * EincMag);
 	sigma_dbsw = 10.*log10(sigma * ko * ko / (4.* Pi * Pi));

//  	fprintf(fsc,
//  		"%5.2f %5.2f %e %e\n"
//  		, theta, phi
//  		, sigma_th_th_dbsw, sigma_ph_th_dbsw );

//   	fprintf(fsc,
//   		"%5.2f %5.2f %e %e\n"
//   		, theta, phi
//   		, sigma_th_th_dbsw, sigma_ph_th_dbsw );

// 	fprintf(fsc,
// 		"%5.2f %5.2f %e %e %e %e\n"
// 		, theta, phi
// 		, (eTheta[i][j]).real(), (eTheta[i][j]).imag()
//                , (ePhi[i][j]).real(),  (ePhi[i][j]).imag());


      if (OPT == 0) {

  	fprintf(fsc,
  		"%5.2f %5.2f %e %e\n"
  		, theta, phi
  		, sigma_th_th_dbsw, sigma_ph_th_dbsw );

      } else {
	fprintf(fsc,
		"%5.2f %5.2f %e %e %e %e\n"
		, theta, phi
		, (eTheta[i][j]).real(), (eTheta[i][j]).imag()
               , (ePhi[i][j]).real(),  (ePhi[i][j]).imag());
      }

};
};
};


  times(&tmsEnd);
  double cpuTime;
  cpuTime = static_cast<double>((tmsEnd.tms_utime + tmsEnd.tms_cutime) -
				(tmsStart.tms_utime + tmsStart.tms_cutime)) /
    static_cast<double>(sysconf(_SC_CLK_TCK));
  int cpuTIME = static_cast<int>(floor(cpuTime + 0.5));
  int hh = cpuTIME / 3600;
  int mm = (cpuTIME - hh * 3600) / 60;
  int ss = (cpuTIME - hh * 3600 - mm * 60);
  cout << "Farfield Time(hh:mm:ss) = " << hh << ":" << mm << ":" << ss << endl;
}

void ReadRegions(char *fname, int **FEM, int **ID, int &numRegions,string** regionName)
{
  char fileName[180];
  char buf[0x200];
  sprintf(fileName, "%s.region", fname);
  ifstream foo;
  foo.open(fileName);
  if (!foo)
    numRegions = 1;
  else 
    foo >> numRegions;
  (*regionName)=new string[numRegions];
  (*FEM) = new int [numRegions];
  (*ID) = new int [numRegions];
  if (!foo) {
    (*ID)[0] = 0;
    (*FEM)[0] = 1;
    return;
  }

  int i, id;
  char solverType[180];
  for (i = 0; i < numRegions; i ++) {
    int id;
    foo >> id;
    foo >> solverType;
    (*ID)[i] = id;
    if (strcmp(solverType, "DDM") == 0) (*FEM)[i] = 1;
    else if (strcmp(solverType, "FEM") == 0) (*FEM)[i] = 1;
    else if (strcmp(solverType, "MOM") == 0) (*FEM)[i] = 0;
    else if (strcmp(solverType, "BEM") == 0) (*FEM)[i] = 0;
    else if (strcmp(solverType, "IE") == 0) (*FEM)[i] = 0;
    else if (strcmp(solverType, "PO") == 0) (*FEM)[i] = 2;
    else if (strcmp(solverType, "FEBI")==0) (*FEM)[i]=3;
    else {
      cout << "unidentified solver type: " << solverType << endl;
      exit(0);
    }
  }
  
  for (i = 0; i < numRegions; i ++) {
	  int id;
	  foo >> id;
	  foo >> buf>> (*regionName)[i];	  
  }
  foo.close();
}