一维层状大地电磁（MT）正演

最近忙先把代码放上吧，图以后再补。杨大神处女作。
C++代码
#include <iostream>
#include <fstream>
#include <time.h>
using namespace std;
const int N = 5; //层数
const int TN = 10; //时间个数
const double PI = 3.0415926535897932384626;
//const double T[100] = { 0.001000,0.001233,0.001520,0.001874,0.002310,0.002848,0.003511,0.004329,0.005337,0.006579,0.008111,0.010000,0.012328,0.015199,0.018738,0.023101,0.028480,0.035112,0.043288,0.053367,0.065793,0.081113,0.100000,0.123285,0.151991,0.187382,0.231013,0.284804,0.351119,0.432876,0.533670,0.657933,0.811131,1.000000,1.232847,1.519911,1.873817,2.310130,2.848036,3.511192,4.328761,5.336699,6.579332,8.111308,10.000000,12.328467,15.199111,18.738174,23.101297,28.480359,35.111917,43.287613,53.366992,65.793322,81.113083,100.000000,123.284674,151.991108,187.381742,231.012970,284.803587,351.119173,432.876128,533.669923,657.933225,811.130831,1000.000000,1232.846739,1519.911083,1873.817423,2310.129700,2848.035868,3511.191734,4328.761281,5336.699231,6579.332247,8111.308308,10000.000000,12328.467394,15199.110830,18738.174229,23101.297001,28480.358684,35111.917342,43287.612811,53366.992312,65793.322466,81113.083079,100000.000000,123284.673944,151991.108295,187381.742286,231012.970008,284803.586844,351119.173422,432876.128108,533669.923121,657933.224658,811130.830790,1000000.000000 };
const double T[10] = { 0.000100,0.001000,0.010000,0.100000,1.000000,10.000000,100.000000,1000.000000,10000.000000,100000.000000 };

//const double rhow[3] = { 50,19*50,19*19*50 };
//const double h[2] = {50,200*50};
const double rhow[5] = { 100,10,1000,1,1000 };
const double h[4] = { 500,500,500,1000 };
class Cfun
{

public:
double L_P(double v, double u, int n);
double L_Q(double v, double u, int n);
double con_Gm(int n, int m);
double R_A(double P, double Q, int n, int m);
double R_B(double P, double Q, int n, int m);
double R_U(double Am, double Bm);
double R_V(double Am, double Bm);
};

double Cfun::L_P(double v, double u,int n)
{
double p1;
double p2;
double p;
p1 = 1.0 - pow(sqrt(rhow[n] / rhow[n - 1]) * u, 2) - pow(sqrt(rhow[n] / rhow[n - 1]) * v, 2);
p2 = pow(1.0 + sqrt(rhow[n] / rhow[n - 1]) * u, 2) + pow(sqrt(rhow[n] / rhow[n - 1]) * v, 2);
p = p1 / p2;
return p;
}
double Cfun::L_Q(double v, double u, int n)
{
double q1;
double q2;
double q;
q1 = -2.0*sqrt(rhow[n] / rhow[n - 1])*v;
q2 = pow(1.0 + sqrt(rhow[n] / rhow[n - 1]) * u, 2) + pow(sqrt(rhow[n] / rhow[n - 1]) * v, 2);
q = q1 / q2;
return q;
}

double Cfun::con_Gm(int n,int m)
{
double g;
g = 4 * PI*h[n] / sqrt(1e7 * rhow[n] * T[m]);
return g;
}

double Cfun::R_A(double P, double Q,int n,int m)
{
double Am;
double Gm = con_Gm(n,m);
Am = exp(-1*Gm)*(P*cos(Gm) - Q*sin(Gm));
return Am;
}
double Cfun::R_B(double P, double Q,int n,int m)
{
double Bm;
double Gm = con_Gm(n,m);
Bm = exp(-1*Gm)*(Q*cos(Gm) + P*sin(Gm));
return Bm;
}
double Cfun::R_U(double Am, double Bm)
{
double u;
u = (1 - Am*Am - Bm*Bm) / (pow(1 + Am, 2) + Bm*Bm);
return u;
}

double Cfun::R_V(double Am, double Bm)
{
double v;
v = (-2 * Bm) / (pow(1 + Am, 2) + Bm*Bm);
return v;
}
int main()
{
time_t start, end;
time(&start);
double p;
double q; 
double Am;
double Bm;
double um;
double vm;
double rt[TN];
//double rabh[TN];
cout << "计算开始……" << endl;

cout << endl;
cout << "**********************************************" << endl;
cout << "* *" << endl;
cout << "* 大地电磁1D_层状正演 *" << endl;
cout << "* Cong Yang *" << endl;
cout << "* CDUT *" << endl;
cout << "* 2017/8/20 *" << endl;
cout << "* *" << endl;
cout << "**********************************************" << endl;
cout << endl;
Cfun Cf;
for (int j = 0; j < TN; j++)
{
um = 1;
vm = 0;
for (int i = N-2; i >= 0; i--)
{
p = Cf.L_P(vm, um, i+1);//m+1
q = Cf.L_Q(vm, um, i+1);//m+1
Am = Cf.R_A(p, q, i, j);//m
Bm = Cf.R_B(p, q, i, j);
um = Cf.R_U(Am, Bm);
vm = Cf.R_V(Am, Bm);
}
rt[j] = rhow[0]*(um*um + vm*vm);
//rabh[j] = sqrt(1e7 * rhow[0] * T[j]) / h[0];
}
ofstream ou("T_rhowT.dat", ios::out);
for (int i = TN-1; i >= 0; i--)
{

//ou << T[i] << "\t" << rabh[i] << endl;
ou << 1/T[i] << "\t" << rt[i] << endl;

}
ou.close();

cout << "结束" << endl;
time(&end);
cout << "用时：" << end - start << "s" << endl;
system("pause");
return 0;
}