二维声波正演 自激自收，二维声波方程有限差分正演，PML吸收边界，

#include <iostream> #include <fstream> #include <cmath> using namespace std; #define M 1000 #define N 120 #define N1 60 #define N2 105 #define X0 500 #define Z0 20 #define Vmax 2600 #define T 1000 int main() { int n, i, j, t; float **u, **w, **p, **q, **r, **pm, **pmu, **pmd, **pm0, **pm1, **pmu1, **pmd1, **pm2, **pmu2, **pmd2, **pm3, **pmu3, **pmd3, **vzu, **vzd, **vz, **ux, **wx, **px, **qx, **rx, **uz, **wz, **pz, **qz, **rz, *dx, *dz; float **ro, **vp, **vs, **c11, **c44, **c13; float dt = 0.001, hx = 5, hz = 5; float c[6] = { 0,1.211243,-0.08972168,0.01384277,-0.00176566,0.0001186795 }; float dxx = log(1 / 0.00001) * 3 * Vmax / (2 * 25 * hx); float dzz = log(1 / 0.00001) * 3 * Vmax / (2 * 25 * hz); ro = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { ro[i] = new float[2 * N + 120]; for (j = 0; j<2 * N1 + 60; j++) ro[i][j] = 1050; for (j = 2 * N1 + 60; j<2 * N2 + 60; j++) ro[i][j] = 1650; for (j = 2 * N2 + 60; j<2 * N + 120; j++) ro[i][j] = 2550; } vp = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { vp[i] = new float[2 * N + 120]; for (j = 0; j<2 * N1 + 60; j++) vp[i][j] = 1500; for (j = 2 * N1 + 60; j<2 * N2 + 60; j++) vp[i][j] = 2000; for (j = 2 * N2 + 60; j<2 * N + 120; j++) vp[i][j] = 2600; } vs = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { vs[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) vs[i][j] = 0; } c11 = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { c11[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) c11[i][j] = vp[i][j] * vp[i][j] * ro[i][j]; } c44 = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { c44[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) c44[i][j] = vs[i][j] * vs[i][j] * ro[i][j]; } c13 = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { c13[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) c13[i][j] = c11[i][j] - 2 * c44[i][j]; } u = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { u[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) u[i][j] = 0; } w = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { w[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) w[i][j] = 0; } p = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { p[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) p[i][j] = 0; } q = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { q[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) q[i][j] = 0; } r = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { r[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) r[i][j] = 0; } pm = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pm[i] = new float[T]; for (j = 0; j<T; j++) pm[i][j] = 0; } pmu = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmu[i] = new float[T]; for (j = 0; j<T; j++) pmu[i][j] = 0; } pmd = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmd[i] = new float[T]; for (j = 0; j<T; j++) pmd[i][j] = 0; } pm0 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pm0[i] = new float[T]; for (j = 0; j<T; j++) pm0[i][j] = 0; } pmu1 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmu1[i] = new float[T]; for (j = 0; j<T; j++) pmu1[i][j] = 0; } pmu2 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmu2[i] = new float[T]; for (j = 0; j<T; j++) pmu2[i][j] = 0; } pmd1 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmd1[i] = new float[T]; for (j = 0; j<T; j++) pmd1[i][j] = 0; } pmd2 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmd2[i] = new float[T]; for (j = 0; j<T; j++) pmd2[i][j] = 0; } pmu3 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmu3[i] = new float[T]; for (j = 0; j<T; j++) pmu3[i][j] = 0; } pmd3 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pmd3[i] = new float[T]; for (j = 0; j<T; j++) pmd3[i][j] = 0; } vzu = new float *[M / 2]; for (i = 0; i<M / 2; i++) { vzu[i] = new float[T]; for (j = 0; j<T; j++) vzu[i][j] = 0; } vzd = new float *[M / 2]; for (i = 0; i<M / 2; i++) { vzd[i] = new float[T]; for (j = 0; j<T; j++) vzd[i][j] = 0; } pm1 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pm1[i] = new float[T]; for (j = 0; j<T; j++) pm1[i][j] = 0; } pm2 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pm2[i] = new float[T]; for (j = 0; j<T; j++) pm2[i][j] = 0; } pm3 = new float *[M / 2]; for (i = 0; i<M / 2; i++) { pm3[i] = new float[T]; for (j = 0; j<T; j++) pm3[i][j] = 0; } vz = new float *[M / 2]; for (i = 0; i<M / 2; i++) { vz[i] = new float[T]; for (j = 0; j<T; j++) vz[i][j] = 0; } ux = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { ux[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) ux[i][j] = 0; } wx = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { wx[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) wx[i][j] = 0; } px = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { px[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) px[i][j] = 0; } qx = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { qx[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) qx[i][j] = 0; } rx = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { rx[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) rx[i][j] = 0; } uz = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { uz[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) uz[i][j] = 0; } wz = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { wz[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) wz[i][j] = 0; } pz = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { pz[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) pz[i][j] = 0; } qz = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { qz[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) qz[i][j] = 0; } rz = new float *[2 * M + 120]; for (i = 0; i<2 * M + 120; i++) { rz[i] = new float[2 * N + 120]; for (j = 0; j<2 * N + 120; j++) rz[i][j] = 0; } dx = new float[2 * M + 120]; dz = new float[2 * N + 120]; for (i = 0; i<2 * M + 120; i++) dx[i] = 0; for (i = 0; i<60; i++) dx[i] = dxx*pow((float(60 - i) / 50), 4); for (i = 2 * M + 60; i<2 * M + 120; i++) dx[i] = dxx*pow((float(i - (2 * M + 59)) / 50), 4); for (j = 0; j<2 * N + 120; j++) dz[j] = 0; for (j = 0; j<60; j++) dz[j] = dzz*pow((float(60 - j) / 50), 4); for (j = 2 * N + 60; j<2 * N + 120; j++) dz[j] = dzz*pow((float(j - (2 * N + 59)) / 50), 4); for (t = 0; t<T; t++) { for (i = 14; i < (M + 30) / 2 - 1; i++) { px[4 * i + 1][2 * Z0 + 60] = px[4 * i + 1][2 * Z0 + 60] + (1 - 2 * pow((3.14 * 35 * (t - 50)*dt), 2))*exp(-pow((3.14 * 35 * (t - 50)*dt), 2)) / 2; pz[4 * i + 1][2 * Z0 + 60] = pz[4 * i + 1][2 * Z0 + 60] + (1 - 2 * pow((3.14 * 35 * (t - 50)*dt), 2))*exp(-pow((3.14 * 35 * (t - 50)*dt), 2)) / 2; p[4 * i + 1][2 * Z0 + 60] = px[4 * i + 1][2 * Z0 + 60] + pz[4 * i + 1][2 * Z0 + 60]; qx[4 * i + 1][2 * Z0 + 60] = qx[4 * i + 1][2 * Z0 + 60] + (1 - 2 * pow((3.14 * 35 * (t - 50)*dt), 2))*exp(-pow((3.14 * 35 * (t - 50)*dt), 2)) / 2; qz[4 * i + 1][2 * Z0 + 60] = qz[4 * i + 1][2 * Z0 + 60] + (1 - 2 * pow((3.14 * 35 * (t - 50)*dt), 2))*exp(-pow((3.14 * 35 * (t - 50)*dt), 2)) / 2; q[4 * i + 1][2 * Z0 + 60] = qx[4 * i + 1][2 * Z0 + 60]