You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
531 lines
15 KiB
531 lines
15 KiB
#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,®ionName);
|
|
|
|
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();
|
|
}
|
|
|