Программа расчета характеристик рассеяния и ослабления электромагнитных волн отдельными сферическими гидрометеорами
-
Просмотров: 571
Программа написана под Builder 6.0 для среды Windows XP.
Текст основного модуля программы.
//------------------------------------------------------------------
//Программа расчета характеристик рассеяния и ослабления //электромагнитных волн отдельными сферическими гидрометеорами
// “Code for electromagnetic scattering by //spherical particles”
/*Входные данные
T1 = n-k*i: комплексный показатель преломления ядра
T2 = n-k*i: комплексный показатель преломления пленки
R=RO - радиус капли и сухого града, см
R=RO+HP если град обводненный, см
HP - толщина пленки, см
LAM - длина волны, см
N, M1 - количество членов ряда в расчетных формулах
N=85; M1=35; если LAM=3.2 см;
N=35; M1=55; если LAM=5.6 - 10.0 см
Выходные данные:
QEXT : поперечник ослабления
QSCA : поперечник полного рассеяния
QRL : поперечник обратного рассеяния */
//------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "Udiggl.h"
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <process.h>
#include <float.h>
#include <complex>
#define len 100
#include <iostream.h>
using namespace std;
//-------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
FILE *AZ;
double SUG (double R, int N, double LAM, int M1);
//сухой град для длин волн от 3,2 до 10 см
//complex <double> T1(1.78,-0.0024);
//complex <double> T2(1.78,-0.0024);
//обводненный град
// LAM=3.2 см
//complex <double> T1(1.78,-0.0024);
//complex <double> T2(7.10,-2.90);
// LAM=5.6 см
complex <double> T1(1.78,-0.0024);
complex <double> T2(8.44,-2.16);
// LAM=10 см
//complex <double> T1(1.78,-0.0024);
//complex <double> T2(9.04,-1.37);
//капли
// LAM=3.2 см
//complex <double> T1(7.10,-2.90);
//complex <double> T2(7.10,-2.90);
// LAM=5.6 см
//complex <double> T1(8.44,-2.16);
//complex <double> T2(8.44,-2.16);
// LAM=10 см
//complex <double> T1(9.04,-1.37);
//complex <double> T2(9.04,-1.37); */
double HP=0.1;
double R,RM;
double LAM;
//-------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
: TForm(Owner)
{
}
//-------------------------------------------------------------------
void __fastcall TForm1::Button1Click(TObject *Sender)
{
int N, M1;
double D;
LAM=5.6;
N=35;
M1=15;
if( (AZ=fopen("az.dat","wt")) == NULL)
{
printf("Error open file AZ.");
exit(1);
};
fprintf (AZ, "внешний радиус,см безразмерный поперечник\n");
R=0.1+HP;
for(int i=0; i<30; i++)
{
D=SUG(R,N,LAM,M1);
fprintf (AZ, "%.5e %.6e\n",R,D);
R=R+0.1;
}
fclose (AZ);
}
double SUG (double R, int N, double LAM, int M1)
{
complex<double> AL1[len],AL2[len],AL3[len],AL4[len];
complex<double> BE1[len],BE2[len],BE3[len],BE4[len];
complex<double> PSI1[len],PSI2[len],PSI3[len],PSI4[len];
complex<double> Y1[len],Y2[len],Y3[len],Y4[len];
complex<double> F1[len],F2[len],F3[len],F4[len];
complex<double> V1[len],V2[len],V3[len],V4[len];
complex<double> KSI1[len],KSI2[len],KSI3[len],KSI4[len];
complex<double> AB[len],BB[len],AM[len],BM[len];
complex<double> Z1,Z2,Z3,T,S2;
double Z4;
double K3,R0;
double DFS,MFS;
int M,L;
double LL1,LL2,LL3,LL4;
double SL,FL,L1,L2,L3,L4,L5,L6,L7,L8;
double X,S1,S3;
double QRL;
M=N-1;
K3=(2.*3.1415926)/LAM;
R0=(R-HP);
T=T1/T2;
Z1=K3*T1*R0;
Z2=K3*T2*R0;
Z3=K3*T2*R;
Z4=K3*R;
F1[N]=0.0;
F1[M]=3.0E-27;
for(L=M; L>=1; L--)
{
L1=L;
F1[L-1]=(2.0*L1+1.0)*F1[L]/Z1-F1[L+1];
}
PSI1[0]=sin(Z1);
for( L=1; L<=M; L++)
{
L2=L;
PSI1[L]=F1[L]*PSI1[0]/F1[0];
Y1[L]=-L2*PSI1[L]/Z1+PSI1[L-1];
}
F2[N]=0.0;
F2[M]=3.0E-27;
for(L=M; L>=1; L--)
{
L3=L;
F2[L-1]=(2.0*L3+1.0)*F2[L]/Z2-F2[L+1];
}
PSI2[0]=sin(Z2);
for( L=1; L<=M; L++)
{
L4=L;
PSI2[L]=F2[L]*PSI2[0]/F2[0];
Y2[L]=-L4*PSI2[L]/Z2+PSI2[L-1];
}
F3[N]=0.0;
F3[M]=3.0E-27;
for(L=M; L>=1; L--)
{
L5=L;
F3[L-1]=(2.0*L5+1.0)*F3[L]/Z3-F3[L+1];
}
PSI3[0]=sin(Z3);
for( L=1; L<=M; L++)
{
L6=L;
PSI3[L]=F3[L]*PSI3[0]/F3[0];
Y3[L]=-L6*PSI3[L]/Z3+PSI3[L-1];
}
F4[N]=0.0;
F4[M]=3.0E-27;
for(L=M; L>=1; L--)
{
L7=L;
F4[L-1]=(2.0*L7+1.0)*F4[L]/Z4-F4[L+1];
}
PSI4[0]=sin(Z4);
for( L=1; L<=M; L++)
{
L8=L;
PSI4[L]=F4[L]*PSI4[0]/F4[0];
Y4[L]=-L8*PSI4[L]/Z4+PSI4[L-1];
}
KSI1[0]=cos(Z1);
KSI1[1]=KSI1[0]/Z1+sin(Z1);
V1[1]=-KSI1[1]/Z1+KSI1[0];
for (L=2; L<=M; L++)
{
LL1=L;
KSI1[L]=(2.0*LL1-1.0)*KSI1[L-1]/Z1-KSI1[L-2];
V1[L]=-LL1*KSI1[L]/Z1+KSI1[L-1];
}
KSI2[0]=cos(Z2);
KSI2[1]=KSI2[0]/Z2+sin(Z2);
V2[1]=-KSI2[1]/Z2+KSI2[0];
for (L=2; L<=M; L++)
{
LL2=L;
KSI2[L]=(2.0*LL2-1.0)*KSI2[L-1]/Z2-KSI2[L-2];
V2[L]=-LL2*KSI2[L]/Z2+KSI2[L-1];
}
KSI3[0]=cos(Z3);
KSI3[1]=KSI3[0]/Z3+sin(Z3);
V3[1]=-KSI3[1]/Z3+KSI3[0];
for (L=2; L<=M; L++)
{
LL3=L;
KSI3[L]=(2.0*LL3-1.0)*KSI3[L-1]/Z3-KSI3[L-2];
V3[L]=-LL3*KSI3[L]/Z3+KSI3[L-1];
}
KSI4[0]=cos(Z4);
KSI4[1]=KSI4[0]/Z4+sin(Z4);
V4[1]=-KSI4[1]/Z4+KSI4[0];
for (L=2; L<=M; L++)
{
LL4=L;
KSI4[L]=(2.0*LL4-1.0)*KSI4[L-1]/Z4-KSI4[L-2];
V4[L]=-LL4*KSI4[L]/Z4+KSI4[L-1];
}
X=K3*R;
S2=0;
complex<double> C(0.0,1.0);
for (L=1; L<=M1; L++)
{
AL1[L]=T*Y2[L]*PSI1[L]-PSI2[L]*Y1[L];
AL2[L]=T*V2[L]*PSI1[L]-KSI2[L]*Y1[L];
BE1[L]=Y2[L]*PSI1[L]/T-PSI2[L]*Y1[L];
BE2[L]=V2[L]*PSI1[L]/T-KSI2[L]*Y1[L];
AL3[L]=AL1[L]*V3[L]-AL2[L]*Y3[L];
AL4[L]=AL1[L]*KSI3[L]-AL2[L]*PSI3[L];
BE3[L]=BE1[L]*V3[L]-BE2[L]*Y3[L];
BE4[L]=BE1[L]*KSI3[L]-BE2[L]*PSI3[L];
AB[L]=(AL3[L]*KSI4[L]-T2*AL4[L]*V4[L])/(AL3[L]*PSI4[L]-T2*AL4[L]*Y4[L]);
BB[L]=(T2*BE3[L]*KSI4[L]-BE4[L]*V4[L])/(T2*BE3[L]*PSI4[L]-BE4[L]*Y4[L]);
AM[L]=1.0/(1.0+C*AB[L]);
BM[L]=1.0/(1.0+C*BB[L]);
double SL;
SL=L;
S1=S1+(2.0*SL+1)*(abs(AM[L])*abs(AM[L])+abs(BM[L])*abs(BM[L]));
S2=S2+pow(-1.0,SL)*(2.0*SL+1.0)*(AM[L]-BM[L]);
S3=S3+(2.0*SL+1)*(real(AM[L]+BM[L]));
}
/*
//поперечник общего рассеяния
QSCA=2.*S1/(X*X);
return (QSCA); */
//поперечник обратного рассеяния
QRL=(abs(S2)*abs(S2))/(X*X);
return (QRL);
/*
//поперечник ослабления
QEXT=2.*S3/(X*X);
return (QEXT); */
}
//-------------------------------------------------------------------
Выходные данные (файл az.dat)
внешний радиус,см безразмерный поперечник
2.00000e-01 7.502621e-03
3.00000e-01 4.065718e-02
4.00000e-01 2.813002e-01
5.00000e-01 8.378652e-01
6.00000e-01 1.540681e+00
7.00000e-01 2.262954e+00
8.00000e-01 2.810482e+00
9.00000e-01 3.009268e+00
1.00000e+00 2.825652e+00
1.10000e+00 2.363339e+00
1.20000e+00 1.771854e+00
1.30000e+00 1.186008e+00
1.40000e+00 6.896357e-01
1.50000e+00 1.606082e-01
1.60000e+00 1.657086e-01
1.70000e+00 6.627160e-01
1.80000e+00 1.223799e+00
1.90000e+00 1.895646e+00
2.00000e+00 2.808344e+00
2.10000e+00 2.981167e+00
2.20000e+00 1.350045e+00
2.30000e+00 7.308839e-01
2.40000e+00 1.730270e-01
2.50000e+00 8.400172e-02
2.60000e+00 5.111780e-01
2.70000e+00 1.822303e+00
2.80000e+00 1.677201e+00
2.90000e+00 1.614138e+00
3.00000e+00 1.501616e+00
3.10000e+00 1.035446e+00
