Skip to content

Commit

Permalink
1: New Interpolation method (use contravarint rather than covarint). …
Browse files Browse the repository at this point in the history 
…2: Modify equatorial wave case. 3: Add the nudging for h.
  • Loading branch information
@Aaron-Hsieh-0129
Aaron-Hsieh-0129 committed Feb 3, 2025
1 parent 47a8851 commit e50aab4
Show file tree
Hide file tree
Showing 9 changed files with 241 additions and 51 deletions.
9 changes: 5 additions & 4 deletions csswm_config.txt
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
CSSWM_OUTPUTPATH=/data/Aaron/TMIF_CSSWM/EquatorialWaveTest
CSSWM_OUTPUTPATH=/data/Aaron/TMIF_CSSWM/Gill/csswm
CSSWM_DT=200
CSSWM_TIMEEND=1200000
CSSWM_OUTPUTSTEP=1
CSSWM_TIMEEND=4320000
CSSWM_OUTPUTSTEP=500
CSSWM_GRAVITY=0.2391
CSSWM_DIFFUSION_KX=200000.
CSSWM_DIFFUSION_KY=200000.
CSSWM_DIFFUSION_TS=0.06
CSSWM_ADDFORCING_TIME=60000
CSSWM_ADDFORCING_TIME=4320000
CSSWM_H_NUDGE_TIME=86400
164 changes: 163 additions & 1 deletion src/bp_wind.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -368,7 +368,7 @@ void CSSWM::BP_wind_interpolation(CSSWM &model) {
#endif
}


/*
void CSSWM::BP_wind_interpolation2(CSSWM &model) {
double B, A1, A2, V1, V2, V3, V4;
Expand Down Expand Up @@ -507,3 +507,165 @@ void CSSWM::BP_wind_interpolation2(CSSWM &model) {
}
#endif
}
*/


void CSSWM::BP_wind_interpolation2(CSSWM &model) {
double B, A1, A2, V1, V2, V3, V4;

int p1, p2, i1, j1, i2, j2, reversed, lonlat;
double uIP, vIP;
int I1, I2_1, I2_2, J1, J2_1, J2_2;
for (int pp = 0; pp < 24; pp++) {
#if defined(SecondOrderSpace)
p1 = match[pp][0], p2 = match[pp][1], i1 = match[pp][2], j1 = match[pp][3], i2 = match[pp][4], j2 = match[pp][5], reversed = match[pp][6], lonlat = match[pp][7];
#elif defined(FourthOrderSpace)
p1 = match_ouTTer[pp][0], p2 = match_ouTTer[pp][1], i1 = match_ouTTer[pp][2], j1 = match_ouTTer[pp][3], i2 = match_ouTTer[pp][4], j2 = match_ouTTer[pp][5], reversed = match_ouTTer[pp][6], lonlat = match_ouTTer[pp][7];
#endif
for (int idx = 0; idx < NX; idx++) {
I1 = i1 == -1 ? idx : i1, J1 = j1 == -1 ? idx : j1;
#if defined(SecondOrderSpace)
I2_1 = i2 == -1 ? reversed ? checkIP[NX-1-idx][0] : checkIP[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP[NY-1-idx][0] : checkIP[idx][0] : j2;
I2_2 = i2 == -1 ? reversed ? checkIP[NX-1-idx][1] : checkIP[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP[NY-1-idx][1] : checkIP[idx][1] : j2;
#elif defined(FourthOrderSpace)
I2_1 = i2 == -1 ? reversed ? checkIP_ouTTer[NX-1-idx][0] : checkIP_ouTTer[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP_ouTTer[NY-1-idx][0] : checkIP_ouTTer[idx][0] : j2;
I2_2 = i2 == -1 ? reversed ? checkIP_ouTTer[NX-1-idx][1] : checkIP_ouTTer[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP_ouTTer[NY-1-idx][1] : checkIP_ouTTer[idx][1] : j2;
#endif
if (lonlat == 0) {
B = model.csswm[p1].lat[I1][J1];
A1 = model.csswm[p2].lat[I2_1][J2_1], A2 = model.csswm[p2].lat[I2_2][J2_2];
// V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2];
// V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2];

V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1];
V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1];
V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2];
V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2];
}
else {
B = model.csswm[p1].lon[I1][J1];
A1 = model.csswm[p2].lon[I2_1][J2_1], A2 = model.csswm[p2].lon[I2_2][J2_2];
// V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2];
// V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2];

V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1];
V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1];
V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2];
V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2];

if (A1 > A2 && (p1 == 0 || p2 == 0)) A2 += 2 * M_PI;
if (A1 > B && B < A2) B += 2 * M_PI;
}

uIP = model.interpolate(A1, A2, V1, V2, B);
vIP = model.interpolate(A1, A2, V3, V4, B);

// TODO:transfer to from contravariant A to covariant B
if (i1 == -1) {
if (j1 == 0) {
#if defined(SecondOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_D[I1][0] * uIP + model.csswm[p1].IP1_D[I1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_D[I1][2] * uIP + model.csswm[p1].IP1_D[I1][3] * vIP;
#elif defined(FourthOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_D[I1][0] * uIP + model.csswm[p1].IP_ouTTer_D[I1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_D[I1][2] * uIP + model.csswm[p1].IP_ouTTer_D[I1][3] * vIP;
#endif
}
else if (j1 == NY-1) {
#if defined(SecondOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_U[I1][0] * uIP + model.csswm[p1].IP1_U[I1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_U[I1][2] * uIP + model.csswm[p1].IP1_U[I1][3] * vIP;
#elif defined(FourthOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_U[I1][0] * uIP + model.csswm[p1].IP_ouTTer_U[I1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_U[I1][2] * uIP + model.csswm[p1].IP_ouTTer_U[I1][3] * vIP;
#endif
}
}
else if (j1 == -1) {
if (i1 == 0) {
#if defined(SecondOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_L[J1][0] * uIP + model.csswm[p1].IP1_L[J1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_L[J1][2] * uIP + model.csswm[p1].IP1_L[J1][3] * vIP;
#elif defined(FourthOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_L[J1][0] * uIP + model.csswm[p1].IP_ouTTer_L[J1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_L[J1][2] * uIP + model.csswm[p1].IP_ouTTer_L[J1][3] * vIP;
#endif
}
else if (i1 == NX-1) {
#if defined(SecondOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_R[J1][0] * uIP + model.csswm[p1].IP1_R[J1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_R[J1][2] * uIP + model.csswm[p1].IP1_R[J1][3] * vIP;
#elif defined(FourthOrderSpace)
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_R[J1][0] * uIP + model.csswm[p1].IP_ouTTer_R[J1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_R[J1][2] * uIP + model.csswm[p1].IP_ouTTer_R[J1][3] * vIP;
#endif
}
}
}
}

#if defined(FourthOrderSpace)
for (int pp = 0; pp < 24; pp++) {
p1 = match_ouTer[pp][0], p2 = match_ouTer[pp][1], i1 = match_ouTer[pp][2], j1 = match_ouTer[pp][3], i2 = match_ouTer[pp][4], j2 = match_ouTer[pp][5], reversed = match_ouTer[pp][6], lonlat = match_ouTer[pp][7];
for (int idx = 0; idx < NX; idx++) {
if (lonlat == 0) {
I1 = i1 == -1 ? idx : i1, J1 = j1 == -1 ? idx : j1;
I2_1 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][0] : checkIP_ouTer[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][0] : checkIP_ouTer[idx][0] : j2;
I2_2 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][1] : checkIP_ouTer[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][1] : checkIP_ouTer[idx][1] : j2;

B = model.csswm[p1].lat[I1][J1];
A1 = model.csswm[p2].lat[I2_1][J2_1], A2 = model.csswm[p2].lat[I2_2][J2_2];
// V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2];
// V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2];

V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1];
V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1];
V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2];
V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2];
}
else {
I1 = i1 == -1 ? idx : i1, J1 = j1 == -1 ? idx : j1;
I2_1 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][0] : checkIP_ouTer[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][0] : checkIP_ouTer[idx][0] : j2;
I2_2 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][1] : checkIP_ouTer[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][1] : checkIP_ouTer[idx][1] : j2;

B = model.csswm[p1].lon[I1][J1];
A1 = model.csswm[p2].lon[I2_1][J2_1], A2 = model.csswm[p2].lon[I2_2][J2_2];
// V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2];
// V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2];

V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1];
V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1];
V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2];
V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2];

if (A1 > A2 && (p1 == 0 || p2 == 0)) A2 += 2 * M_PI;
if (A1 > B && B < A2) B += 2 * M_PI;
}

uIP = model.interpolate(A1, A2, V1, V2, B);
vIP = model.interpolate(A1, A2, V3, V4, B);

if (i1 == -1) {
if (j1 == 1) {
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_D[I1][0] * uIP + model.csswm[p1].IP_ouTer_D[I1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_D[I1][2] * uIP + model.csswm[p1].IP_ouTer_D[I1][3] * vIP;
}
else if (j1 == NY-2) {
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_U[I1][0] * uIP + model.csswm[p1].IP_ouTer_U[I1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_U[I1][2] * uIP + model.csswm[p1].IP_ouTer_U[I1][3] * vIP;
}
}
else if (j1 == -1) {
if (i1 == 1) {
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_L[J1][0] * uIP + model.csswm[p1].IP_ouTer_L[J1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_L[J1][2] * uIP + model.csswm[p1].IP_ouTer_L[J1][3] * vIP;
}
else if (i1 == NX-2) {
model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_R[J1][0] * uIP + model.csswm[p1].IP_ouTer_R[J1][1] * vIP;
model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_R[J1][2] * uIP + model.csswm[p1].IP_ouTer_R[J1][3] * vIP;
}
}
}
}
#endif
}
4 changes: 4 additions & 0 deletions src/construction.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ class CSSWM {
double du[NX][NY][2];
double dv[NX][NY][2];
#endif

double f[NX][NY];
};

CSSWM();
Expand Down Expand Up @@ -79,6 +81,7 @@ class CSSWM {
double diffusion_ky = KY;
double diffusion_ts = TIMETS;
double addforcingtime = ADDFORCINGTIME;
double csswm_h_nudge_time = CSSWM_H_NUDGE_TIME;


// ***********************************************************************************
Expand Down Expand Up @@ -195,6 +198,7 @@ class CSSWM {
public:
static void DiffusionAll(CSSWM &);
static void timeFilterAll(CSSWM &);
static void NudgeH(CSSWM &);
};


Expand Down
2 changes: 2 additions & 0 deletions src/define.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#define ALPHA0 (0)
// #define Advection
// #define GravityWave
// #define DeformationalFlow
// #define SteadyGeostrophy
// #define Barotropic
// #define Mountain
Expand All @@ -49,5 +50,6 @@
#define TIMETS (0.06)

#define ADDFORCINGTIME (100. * 60.)
#define CSSWM_H_NUDGE_TIME (0.)

#endif
20 changes: 16 additions & 4 deletions src/initialField.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@ void CSSWM::Init::Init2d(CSSWM & model) {

#ifdef DeformationalFlow
model.csswm[p].hp[i][j] = DeformationalFlowH(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]);
model.csswm[p].up[i][j] = 0.;
model.csswm[p].vp[i][j] = 0.;
#endif

#ifdef GravityWave
Expand All @@ -48,6 +50,8 @@ void CSSWM::Init::Init2d(CSSWM & model) {
(model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * SteadyGeostrophyV(model.csswm[p].lon[i][j]);
model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * SteadyGeostrophyU(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]) +
(model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * SteadyGeostrophyV(model.csswm[p].lon[i][j]);

model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0));
#endif

#ifdef Barotropic
Expand All @@ -57,6 +61,7 @@ void CSSWM::Init::Init2d(CSSWM & model) {
(model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * 0;
model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * BarotropicU(model.csswm[p].lat[i][j]) +
(model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * 0;
model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0));
#endif

#ifdef Mountain
Expand All @@ -65,6 +70,7 @@ void CSSWM::Init::Init2d(CSSWM & model) {
(model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * MountainV(model.csswm[p].lon[i][j]);
model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * MountainU(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]) +
(model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * MountainV(model.csswm[p].lon[i][j]);
model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0));
#endif

#ifdef RossbyHaurwitz
Expand All @@ -73,10 +79,11 @@ void CSSWM::Init::Init2d(CSSWM & model) {
(model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * RossbyHaurwitzV(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]);
model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * RossbyHaurwitzU(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]) +
(model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * RossbyHaurwitzV(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]);
model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0));
#endif

#ifdef EquatorialWave
if (p == 0) {
if (p == 1) {
model.csswm[p].h_forcing[i][j] = EquatorialWaveH(model.csswm[p].x[i][j], model.csswm[p].y[i][j]);
}
else {
Expand All @@ -86,6 +93,11 @@ void CSSWM::Init::Init2d(CSSWM & model) {
model.csswm[p].hp[i][j] = 10454.608791605699 + model.csswm[p].h_forcing[i][j];
model.csswm[p].up[i][j] = 0.;
model.csswm[p].vp[i][j] = 0.;

// double f0 = 0;
// double beta = 2.5 * 10E-11;
// model.csswm[p].f[i][j] = f0 + beta * model.csswm[p].lat[i][j] * 180. / M_PI * (111000.) / 11.5;
model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0));
#endif

#if defined(Uniform)
Expand Down Expand Up @@ -302,7 +314,7 @@ double CSSWM::Init::RossbyHaurwitzV(double lon, double lat) {
}

double CSSWM::Init::EquatorialWaveH(double x, double y) {
double h0 = 0.08;
double a = 100000., b = 100000.;
return -h0 * exp(-(x*x / (2*(a*a)) + y*y / (2*(b*b))));
double h0 = 5E-2;
double a = 400. * 1E3, b = 200. * 1E3;
return -h0 * exp(-(x*x / ((a*a)) + y*y / ((b*b))));
}
Loading

0 comments on commit e50aab4

Please sign in to comment.