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.
150 lines
2.9 KiB
150 lines
2.9 KiB
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <iostream>
|
|
#include "integral.h"
|
|
#include "transform.h"
|
|
#include "global.h"
|
|
|
|
cVtr Near2FarTransform(int numAbcTri,
|
|
double k,
|
|
vtr O,
|
|
double theta,
|
|
double phi,
|
|
regcvtr *regJ,
|
|
regcvtr *regM,
|
|
tri3d **triArray)
|
|
{
|
|
int i;
|
|
cVtr E, A, F;
|
|
vtr rvtr;
|
|
|
|
E.reset();
|
|
A.reset(); F.reset();
|
|
rvtr = RDirection(theta, phi);
|
|
for (i = 0; i < numAbcTri; i ++) {
|
|
tri3d *fc;
|
|
cVtr triA, triF;
|
|
|
|
fc = triArray[i];
|
|
TriGaussIntegral(fc, k, rvtr, O, regJ[i], regM[i], 7, &triA, &triF);
|
|
|
|
A = A + triA;
|
|
F = F + triF;
|
|
}
|
|
E = A + F * rvtr;
|
|
E = E * k;
|
|
|
|
return E;
|
|
}
|
|
|
|
cVtr Near2FarTransform(int numAbcTri,
|
|
double k,
|
|
vtr O,
|
|
double theta,
|
|
double phi,
|
|
regcvtr *regJ,
|
|
tri3d **triArray)
|
|
{
|
|
int i;
|
|
cVtr E, A, F;
|
|
vtr rvtr;
|
|
|
|
A.reset(); F.reset();
|
|
rvtr = RDirection(theta, phi);
|
|
for (i = 0; i < numAbcTri; i ++) {
|
|
tri3d *fc;
|
|
cVtr triA, triF;
|
|
|
|
fc = triArray[i];
|
|
TriGaussIntegral(fc, k, rvtr, O, regJ[i], 7, &triA);
|
|
|
|
A = A + triA;
|
|
}
|
|
|
|
E = A;
|
|
E = E * k;
|
|
|
|
return E;
|
|
}
|
|
|
|
void TriGaussIntegral(tri3d *fc,
|
|
double k,
|
|
vtr rvtr,
|
|
vtr O,
|
|
regcvtr regJ,
|
|
regcvtr regM,
|
|
int np,
|
|
cVtr *J,
|
|
cVtr *M)
|
|
{
|
|
int i;
|
|
double zeta0, zeta1, zeta2, weight;
|
|
cVtr cvJ, cvM;
|
|
|
|
J->reset();
|
|
M->reset();
|
|
for (i = 0; i < np; i ++) {
|
|
vtr rp;
|
|
double x, y, z;
|
|
Complex cval;
|
|
|
|
GetFormula(np, i, &zeta0, &zeta1, &zeta2, &weight);
|
|
GetCartesian(fc, zeta0, zeta1, zeta2, &x, &y, &z);
|
|
rp.setvtr(x, y, z);
|
|
rp = rp - O;
|
|
|
|
double val = k * dotP(rp, rvtr);
|
|
cval = Complex (cos(val), sin(val));
|
|
|
|
cvJ = GetField(regJ, zeta0, zeta1, zeta2);
|
|
cvM = GetField(regM, zeta0, zeta1, zeta2);
|
|
|
|
(*J) = (*J) + (cvJ * cval) * (weight * No);
|
|
(*M) = (*M) + (cvM * cval) * weight;
|
|
}
|
|
|
|
double area;
|
|
// area = fc->GETAREA();
|
|
area = Area(fc);
|
|
(*J) = (*J) * fabs(area);
|
|
(*M) = (*M) * fabs(area);
|
|
}
|
|
|
|
void TriGaussIntegral(tri3d *fc,
|
|
double k,
|
|
vtr rvtr,
|
|
vtr O,
|
|
regcvtr regJ,
|
|
int np,
|
|
cVtr *J)
|
|
{
|
|
int i;
|
|
double zeta0, zeta1, zeta2, weight;
|
|
cVtr cvJ;
|
|
|
|
J->reset();
|
|
for (i = 0; i < np; i ++) {
|
|
vtr rp;
|
|
double x, y, z;
|
|
Complex cval;
|
|
|
|
GetFormula(np, i, &zeta0, &zeta1, &zeta2, &weight);
|
|
GetCartesian(fc, zeta0, zeta1, zeta2, &x, &y, &z);
|
|
rp.setvtr(x, y, z);
|
|
rp = rp - O;
|
|
|
|
double val = k * dotP(rp, rvtr);
|
|
cval = Complex (cos(val), sin(val));
|
|
|
|
cvJ = GetField(regJ, zeta0, zeta1, zeta2);
|
|
|
|
(*J) = (*J) + (cvJ * cval) * (weight * No);
|
|
}
|
|
|
|
double area;
|
|
area = Area(fc);
|
|
(*J) = (*J) * fabs(area);
|
|
}
|
|
|
|
|
|
|
|
|