📄 total_power3.c
字号:
if( (flagP2 == 2))
{
ind2 = ind*6;
epsreal = material[4 + ind2];
epsimag = -60.0*lambda*material[5 + ind2];
tmp = 1.0/dj;
Vd = planes[12 + ind1]*P1P2x*tmp + uy*P1P2y*tmp + planes[14 + ind1]*P1P2z*tmp;
vs = fabs(Vd);
vs2 = vs*vs;
x1 = epsreal - (1.0 - vs2);
y1 = epsimag;
r1 = sqrt(x1*x1 + y1*y1);
x2 = sqrt((r1 + x1)*0.5);
y2 = y1/(2.0*x2);
if(fabs(uy) > 0.5)
{
x3 = epsreal*vs - x2;
y3 = epsimag*vs - y2;
x4 = epsreal*vs + x2;
y4 = epsimag*vs + y2;
r2 = 1.0/(x4*x4 + y4*y4);
x5 = (x3*x4 + y3*y4)*r2;
y5 = (y3*x4 - x3*y4)*r2;
loss_jreal = sqrt(loss_factor*(1.0 - (x5*x5 + y5*y5)));
}
else
{
x3 = vs - x2;
y3 = -y2;
x4 = vs + x2;
y4 = y2;
r2 = 1.0/(x4*x4 + y4*y4);
x5 = (x3*x4 + y3*y4)*r2;
y5 = (y3*x4 - x3*y4)*r2;
loss_jreal = sqrt(loss_factor*(1.0 - (x5*x5 + y5*y5)));
}
}
distance[currentpath] += dj;
tmp1 = (lossreal*loss_jreal - lossimag*loss_jimag);
tmp2 = (lossreal*loss_jimag + lossimag*loss_jreal);
lossreal = tmp1;
lossimag = tmp2;
}
lossfac[currentpath] = (Ao*sqrt(lossreal*lossreal + lossimag*lossimag))/(distance[currentpath]*cte);
}
/* ----------------------------------------------------------------------------- */
void fct_calimage(double *TX0 , double *planes , int i , double *TX_i)
{
int id = i*22;
register double t_imi , ux = planes[12 + id] , uy = planes[13 + id] , uz = planes[14 + id];
register double TX0x = TX0[0] , TX0y = TX0[1] , TX0z = TX0[2];
t_imi = -2.0*(TX0x*ux + TX0y*uy + TX0z*uz + planes[15 + id]);
TX_i[0] = TX0x + ux*t_imi;
TX_i[1] = TX0y + uy*t_imi;
TX_i[2] = TX0z + uz*t_imi;
}
/* ----------------------------------------------------------------------------- */
bool fct_calrpt(double *P1, double *P2, double *planes , int i , double *Phit)
{
int id = i*22;
double ux = planes[12 + id], uy = planes[13 + id] , uz = planes[14 + id];
register double P1x = P1[0], P1y = P1[1] , P1z = P1[2];
double ti_max, ti_maxtol , Vdi , P1P2x , P1P2y , P1P2z , Rdx , Rdy , Rdz , ti , cte;
P1P2x = P2[0] - P1x;
P1P2y = P2[1] - P1y;
P1P2z = P2[2] - P1z;
ti_max = sqrt(P1P2x*P1P2x + P1P2y*P1P2y + P1P2z*P1P2z);
ti_maxtol = ti_max - t_tol;
if (ti_max > t_tol)
{
cte = 1.0/ti_max;
Rdx = P1P2x*cte;
Rdy = P1P2y*cte;
Rdz = P1P2z*cte;
Vdi = ux*Rdx + uy*Rdy + uz*Rdz;
if (Vdi != 0.0)
{
ti = - (ux*P1x + uy*P1y + uz*P1z + planes[15 + id]) / Vdi;
if ((t_tol < ti) && (ti < ti_maxtol))
{
Phit[0] = P1x + Rdx*ti;
Phit[1] = P1y + Rdy*ti;
Phit[2] = P1z + Rdz*ti;
if ((planes[16 + id] <= Phit[0]) && (Phit[0] <= planes[17 + id]) && (planes[18 + id] <= Phit[1]) && (Phit[1]<=planes[19 + id]) && (planes[20 + id] <= Phit[2]) && (Phit[2] <= planes[21 + id]))
{
return true;
}
}
}
}
return false;
}
/* ----------------------------------------------------------------------------- */
double * fct_caltpts(double *P1, double *P2, double *planes, int nplanes , double *material , int *nthits , double *Phit_nt)
{
int number = 0 , nt , ntd , jd , i , ind , id;
double tmax_nt, tmax_ntol;
double *result=NULL;
double Vd_nt , P1P2x , P1P2y , P1P2z , Rdx_nt , Rdy_nt , Rdz_nt , t_nt , tmp ;
double P1x = P1[0], P1y = P1[1] , P1z = P1[2];
double ux_nt , uy_nt, uz_nt;
double *tab_pts=NULL , *sorted_tab_pts=NULL;
int *index=NULL;
bool already_there = false;
P1P2x = P2[0] - P1x;
P1P2y = P2[1] - P1y;
P1P2z = P2[2] - P1z;
tmax_nt = sqrt(P1P2x*P1P2x + P1P2y*P1P2y + P1P2z*P1P2z);
tmax_ntol = tmax_nt - t_tol;
if (tmax_nt > t_tol)
{
tmp = 1.0/tmax_nt ;
Rdx_nt = P1P2x*tmp;
Rdy_nt = P1P2y*tmp;
Rdz_nt = P1P2z*tmp;
for (nt = 0 ; nt < nplanes ; nt++)
{
ntd = nt*22;
ux_nt = planes[12 + ntd];
uy_nt = planes[13 + ntd];
uz_nt = planes[14 + ntd];
Vd_nt = ux_nt*Rdx_nt + uy_nt*Rdy_nt + uz_nt*Rdz_nt;
if (Vd_nt != 0.0)
{
t_nt = - (ux_nt*P1x + uy_nt*P1y + uz_nt*P1z + planes[15 + ntd]) / Vd_nt;
if ((t_tol < t_nt) && (t_nt < tmax_ntol))
{
Phit_nt[0] = P1x + Rdx_nt*t_nt;
Phit_nt[1] = P1y + Rdy_nt*t_nt;
Phit_nt[2] = P1z + Rdz_nt*t_nt;
if ((planes[16 + ntd] <= Phit_nt[0]) && (Phit_nt[0] <= planes[17 + ntd]) && (planes[18 + ntd] <= Phit_nt[1]) && (Phit_nt[1] <= planes[19 + ntd]) && ( planes[20 + ntd] <= Phit_nt[2]) && (Phit_nt[2] <= planes[21 + ntd]))
{
if(number == 0)
{
jd = number*6;
number++;
tab_pts = (double *)realloc((double *)tab_pts , (jd + 6)*sizeof(double));
tab_pts[0 + jd] = Phit_nt[0];
tab_pts[1 + jd] = Phit_nt[1];
tab_pts[2 + jd] = Phit_nt[2];
tab_pts[3 + jd] = nt;
tab_pts[4 + jd] = 2.0;
tab_pts[5 + jd] = t_nt;
}
else
{
already_there = false;
for (i = 0 ; i < number ; i++)
{
id = 6*i;
if( (fabs(tab_pts[0 + id] - Phit_nt[0]) < t_tol) && (fabs(tab_pts[1 + id] - Phit_nt[1]) < t_tol) && (fabs(tab_pts[2 + id] - Phit_nt[2]) < t_tol) )
{
already_there = true;
}
}
if(already_there == false)
{
jd = number*6;
number++;
tab_pts = (double *)realloc((double *)tab_pts , (jd + 6)*sizeof(double));
tab_pts[0 + jd] = Phit_nt[0];
tab_pts[1 + jd] = Phit_nt[1];
tab_pts[2 + jd] = Phit_nt[2];
tab_pts[3 + jd] = nt;
tab_pts[4 + jd] = 2.0;
tab_pts[5 + jd] = t_nt;
}
}
}
}
}
}
if(number > 0)
{
sorted_tab_pts = (double *)malloc(number*sizeof(double));
index = (int *)malloc(number*sizeof(int));
for (i = 0 ; i < number ; i++)
{
index[i] = i;
sorted_tab_pts[i] = tab_pts[5 + i*6];
}
qsindex( sorted_tab_pts, index , 0 , number - 1 );
result = (double *)malloc(number*5*sizeof(double));
for (i = 0 ; i < number ; i++)
{
ind = 6*index[i];
id = 5*i;
result[0 + id] = tab_pts[0 + ind];
result[1 + id] = tab_pts[1 + ind];
result[2 + id] = tab_pts[2 + ind];
result[3 + id] = tab_pts[3 + ind];
result[4 + id] = tab_pts[4 + ind];
}
}
}
free(tab_pts);
free(sorted_tab_pts);
free(index);
nthits[0] = number;
return result;
}
/* ----------------------------------------------------------------------------- */
void qsindex (double *a, int *index , int lo, int hi)
{
// lo is the lower index, hi is the upper index
// of the region of array a that is to be sorted
int i=lo, j=hi , ind;
double x=a[(lo+hi)/2] , h;
// partition
do
{
while (a[i]<x) i++;
while (a[j]>x) j--;
if (i<=j)
{
h = a[i];
a[i] = a[j];
a[j] = h;
ind = index[i];
index[i] = index[j];
index[j] = ind;
i++;
j--;
}
}
while (i<=j);
// recursion
if (lo<j) qsindex(a , index , lo , j);
if (i<hi) qsindex(a , index , i , hi);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -