最近忙先把代码放上吧,图以后再补。杨大神处女作。 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; }
文章评论