00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <float.h>
00026 #include <math.h>
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include "defs.h"
00030 #include "athena.h"
00031 #include "globals.h"
00032 #include "prototypes.h"
00033
00034
00035 static const Real mbar = (1.27)*(1.6733e-24);
00036 static const Real kb = 1.380658e-16;
00037 static const Real pi = 3.14159265358979323846;
00038
00039
00040
00041
00042
00043
00044 static double ran2(long int *idum);
00045 static Real logd(const Grid *pG, const int i, const int j, const int k);
00046
00047
00048
00049
00050
00051 void problem(Grid *pGrid, Domain *pDomain)
00052 {
00053 int i=0,j=0,k=0,n=0,m=0,l=0;
00054 int is,ie,js,je,ks,ke,iprob,
00055 nkx=16,nky=8,nkz=8,nx1,nx2,nx3;
00056 Real P_k, n0, T0, kappa,chi, x1L, x2L, x3L, x12L, num,
00057 drho, dp,
00058 krho, krho_i, krho_ij, angle, phi, phi2, phi3, ksi, tmp,
00059 rho0, rho1, n_anal, r, r_ij, r_proj,
00060 beta, b0, Btot=0.0,Bsum=0.0,
00061 ***az,***ax,***ay,amp,kx[16],ky[8],kz[8],xa,ya,za;
00062 long int iseed = -1, rseed;
00063
00064 is = pGrid->is; ie = pGrid->ie;
00065 js = pGrid->js; je = pGrid->je;
00066 ks = pGrid->ks; ke = pGrid->ke;
00067
00068
00069 P_k = par_getd("problem","P_k");
00070 n0 = par_getd("problem","n0");
00071 x1L = par_getd("grid","x1max");
00072 x2L = par_getd("grid","x2max");
00073 x3L = par_getd("grid","x3max");
00074 iprob = par_geti("problem","iprob");
00075
00076 #ifdef MHD
00077 beta = par_getd("problem","beta");
00078
00079 b0 = sqrt(2.0*P_k*kb/beta);
00080 #endif
00081 #ifdef ISOTROPIC_CONDUCTION
00082 kappa = par_getd("problem","kappa");
00083 kappa_T = (mbar/kb)*kappa;
00084 #endif
00085 #ifdef ANISOTROPIC_CONDUCTION
00086 chi = par_getd("problem","chi");
00087 chi_C = (mbar/kb)*chi;
00088 #endif
00089
00090
00091
00092 T0 = P_k / n0;
00093 rho0 = n0 * mbar;
00094
00095
00096
00097
00098
00099 if(iprob==1) {
00100 drho = par_getd("problem","drho");
00101 n_anal = par_getd("problem","n_anal");
00102 num = par_getd("problem","num");
00103
00104 krho_i = 2.0 * pi * num / (float)(ie-is+1);
00105 krho = 2.0 * pi * num / x1L;
00106 rho1 = rho0 * drho;
00107 for (k=ks; k<=ke; k++) {
00108 for (j=js; j<=je; j++) {
00109 for (i=is; i<=ie; i++) {
00110 pGrid->U[k][j][i].d = rho0 + rho1*cos(krho_i*(float)(i-is));
00111 pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d * n_anal * rho1 * sin(krho_i*(float)(i-is)) / rho0 / krho * (-1.0);
00112 pGrid->U[k][j][i].M2 = 0.0;
00113 pGrid->U[k][j][i].M3 = 0.0;
00114 pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 - n_anal*n_anal/ (krho * krho * Gamma_1) * rho1*cos(krho_i*(float)(i-is))
00115 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00116 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00117 }
00118 }
00119 }
00120 }
00121
00122
00123
00124
00125
00126
00127 if(iprob==2) {
00128 angle = atan(x2L/x1L);
00129 x12L = pow(x1L*x1L + x2L*x2L, 0.5);
00130 r_ij = pow(pow((float)(ie-is),2)+pow((float)(je-js),2),0.5);
00131 krho = 2.0 * pi * num / x12L;
00132 krho_ij = 2.0 * pi * num / r_ij;
00133 rho1 = rho0 * drho;
00134 for (k=ks; k<=ke; k++) {
00135 for (j=js; j<=je; j++) {
00136 for (i=is; i<=ie; i++) {
00137 r = pow(pow((float)(i-is),2)+pow((float)(j-js),2),0.5);
00138 if ((i-is) == 0) {
00139 phi = 0;
00140 }
00141 else{
00142 phi = atan((float)(j-js)/(float)(i-is));
00143 }
00144 r_proj = r * sin(phi + angle) / sin ( 2 * angle);
00145 pGrid->U[k][j][i].d = rho0 + rho1 * cos(krho_ij * r_proj);
00146 pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d * n_anal * rho1 * sin(krho_ij * r_proj) / rho0 / krho * cos(angle) * (-1.0);
00147 pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d * n_anal * rho1 * sin(krho_ij * r_proj) / rho0 / krho * sin(angle) * (-1.0);
00148 pGrid->U[k][j][i].M3 = 0.0;
00149 pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1
00150 - n_anal*n_anal/ (krho * krho * Gamma_1) * rho1 * cos(krho_ij * r_proj)
00151 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00152 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00153 }
00154 }
00155 }
00156 }
00157
00158
00159
00160
00161
00162
00163 if(iprob==3) {
00164 dp = par_getd("problem","dp");
00165 angle = atan(x2L/x1L);
00166 for (k=ks; k<=ke; k++) {
00167 for (j=js; j<=je; j++) {
00168 for (i=is; i<=ie; i++) {
00169 pGrid->U[k][j][i].d = rho0;
00170 pGrid->U[k][j][i].M1 = 0.0;
00171 pGrid->U[k][j][i].M2 = 0.0;
00172 pGrid->U[k][j][i].M3 = 0.0;
00173 pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00174 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00175 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00176 #ifdef MHD
00177 pGrid->B1i[k][j][i] = b0 * cos(angle);
00178 if (i == ie) pGrid->B1i[k][j][i+1] = b0 * cos(angle);
00179 pGrid->U[k][j][i].B1c = b0 * cos(angle);
00180
00181 pGrid->B2i[k][j][i] = b0 * sin(angle);
00182 if (j==je) pGrid->B2i[k][j+1][i] = b0 * sin(angle);
00183 pGrid->U[k][j][i].B2c = b0 * sin(angle);
00184
00185 pGrid->B3i[k][j][i] = 0.0;
00186 if (k==ke && pGrid->Nx3 > 1) pGrid->B3i[k+1][j][i] = 0.0;
00187 pGrid->U[k][j][i].B3c = 0.0;
00188 pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00189 + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00190 #endif
00191 }
00192 }
00193 }
00194 }
00195
00196
00197
00198
00199
00200 if(iprob==4){
00201
00202 nx1 = (ie-is)+1 + 2*nghost;
00203 nx2 = (je-js)+1 + 2*nghost;
00204 if ((nx1 == 1) || (nx2 == 1)) {
00205 ath_error("[TI-tangledB]: This problem can only be run with Nx1>1\n");
00206 }
00207 if ((az = (Real***)calloc_3d_array(nx3,nx2, nx1, sizeof(Real))) == NULL) {
00208 ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00209 }
00210
00211
00212 for (n=0; n<=nkx-1; n++){
00213 kx[n] = (1.+(float)(n))*(2.0*pi)/2.0;
00214 }
00215 kx[nkx-1] = kx[0] * (-1.0);
00216 for (m=0; m<=nky-1; m++){
00217 ky[m] = (1.+(float)(m))*(2.0*pi);
00218 }
00219 ky[nky-1] = ky[0] * (-1.0);
00220
00221 for (k=ks; k=ks; k++) {
00222 for (j=js; j<=je+1; j++) {
00223 for (i=is; i<=ie+1; i++) {
00224 az[k][j][i] = 0.0;
00225 }
00226 }
00227 }
00228 for (n=0; n<=nkx-1; n++){
00229 for (m=0; m<=nky-1; m++){
00230 phi = 2.0 * pi * ran2(&rseed);
00231
00232 amp = (1.0/sqrt( kx[n]*kx[n] + ky[m]*ky[m]));
00233 for (j=js; j<=je+1; j++) {
00234 ya=(float)(j-js)*(1.0/(float)(je-js+1));
00235
00236
00237
00238 for (i=is; i<=ie+1; i++) {
00239 xa=(float)(i-is)*(2.0/(float)(ie-is+1));
00240 az[ks][j][i] += amp*sin(kx[n]*xa + ky[m]*ya + phi);
00241 }
00242 }
00243 }
00244 }
00245
00246
00247 #ifdef MHD
00248 dp = par_getd("problem","dp");
00249 for (k=ks; k<=ke; k++) {
00250 for (j=js; j<=je; j++) {
00251 for (i=is; i<=ie; i++) {
00252 pGrid->U[k][j][i].d = rho0;
00253 pGrid->U[k][j][i].M1 = 0.0;
00254 pGrid->U[k][j][i].M2 = 0.0;
00255 pGrid->U[k][j][i].M3 = 0.0;
00256 pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00257 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00258 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00259 pGrid->B1i[k][j][i] = (az[ks][j+1][i] - az[ks][j][i]);
00260 pGrid->B2i[k][j][i] =-(az[ks][j][i+1] - az[ks][j][i]);
00261
00262 pGrid->B3i[k][j][i] = 0.0;
00263 if (k==ke && pGrid->Nx3 > 1) pGrid->B3i[k+1][j][i] = 0.0;
00264
00265 Bsum += (pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00266 pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00267 }
00268 }
00269 }
00270 Btot = sqrt(Bsum/((float)(je-js+1) * (float)(ie-is+1)));
00271 printf("Btot=%f\n",Btot);
00272 Bsum = 0.0;
00273
00274 for (k=ks; k<=ke; k++) {
00275 for (j=js; j<=je; j++) {
00276 for (i=is; i<=ie; i++) {
00277 pGrid->B1i[k][j][i] = pGrid->B1i[k][j][i] * b0 / Btot;
00278 if (i == ie) pGrid->B1i[k][j][i+1] = pGrid->B1i[k][j][is];
00279 pGrid->B2i[k][j][i] = pGrid->B2i[k][j][i] * b0 / Btot;
00280 if (j == je) pGrid->B2i[k][j+1][i] = pGrid->B2i[k][js][i];
00281 Bsum += (pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00282 pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00283 }
00284 }
00285 }
00286
00287 for (k=ks; k<=ke; k++) {
00288 for (j=js; j<=je; j++) {
00289 for (i=is; i<=ie; i++) {
00290 pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00291 pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00292 pGrid->U[k][j][i].B3c = 0.0;
00293 pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00294 + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00295 }
00296 }
00297 }
00298 Btot = Bsum/((float)(je-js+1) * (float)(ie-is+1));
00299 printf("B^2tot,b0^2=%e %e \n",Btot,b0*b0);
00300 #endif //MHD
00301 }
00302
00303
00304
00305
00306
00307
00308 if(iprob==5) {
00309 dp = par_getd("problem","dp");
00310 angle = atan(x2L/x1L);
00311 x12L = pow(x1L*x1L + x2L*x2L, 0.5);
00312 ksi = atan(x3L/x12L);
00313 for (k=ks; k<=ke; k++) {
00314 for (j=js; j<=je; j++) {
00315 for (i=is; i<=ie; i++) {
00316 pGrid->U[k][j][i].d = rho0;
00317 pGrid->U[k][j][i].M1 = 0.0;
00318 pGrid->U[k][j][i].M2 = 0.0;
00319 pGrid->U[k][j][i].M3 = 0.0;
00320 pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00321 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00322 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00323 #ifdef MHD
00324 pGrid->B1i[k][j][i] = b0 * cos(angle) * cos(ksi);
00325 if (i == ie) pGrid->B1i[k][j][i+1] = b0 * cos(angle) * cos(ksi);
00326 pGrid->U[k][j][i].B1c = b0 * cos(angle) * cos(ksi);
00327
00328 pGrid->B2i[k][j][i] = b0 * sin(angle) * cos(ksi);
00329 if (j==je) pGrid->B2i[k][j+1][i] = b0 * sin(angle) * cos(ksi);
00330 pGrid->U[k][j][i].B2c = b0 * sin(angle) * cos(ksi);
00331
00332 pGrid->B3i[k][j][i] = b0 * sin(ksi);
00333 if (k==ke && pGrid->Nx3 > 1) pGrid->B3i[k+1][j][i] = b0 * sin(ksi);
00334 pGrid->U[k][j][i].B3c = b0 * sin(ksi);
00335 pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00336 + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00337 #endif
00338 }
00339 }
00340 }
00341 }
00342
00343
00344
00345
00346
00347
00348 if(iprob==6){
00349 nx1 = (ie-is)+1 + 2*nghost;
00350 nx2 = (je-js)+1 + 2*nghost;
00351 nx3 = (ke-ks)+1 + 2*nghost;
00352 if ((nx1 == 1) || (nx2 == 1) || (nx3 == 1)) {
00353 ath_error("[TI-tangledB]: This problem can only be run with Nx1>1\n");
00354 }
00355 if ((ay = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00356 ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00357 }
00358 if ((az = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00359 ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00360 }
00361 if ((ax = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00362 ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00363 }
00364
00365
00366 for (n=0; n<=nkx-1; n++){
00367 kx[n] = (1.+(float)(n))*(2.0*pi)/2.0;
00368 }
00369 kx[nkx-1] = kx[0] * (-1.0);
00370 for (m=0; m<=nky-1; m++){
00371 ky[m] = (1.+(float)(m))*(2.0*pi);
00372 }
00373 ky[nky-1] = ky[0] * (-1.0);
00374 for (l=0; l<=nkz-1; l++){
00375 kz[l] = (1.+(float)(l))*(2.0*pi);
00376 }
00377 kz[nkz-1] = kz[0] * (-1.0);
00378
00379 for (k=ks; k<=ke+1; k++) {
00380 for (j=js; j<=je+1; j++) {
00381 for (i=is; i<=ie+1; i++) {
00382 ax[k][j][i] = 0.0;
00383 ay[k][j][i] = 0.0;
00384 az[k][j][i] = 0.0;
00385 }
00386 }
00387 }
00388
00389 for (n=0; n<=nkx-1; n++){
00390 for (m=0; m<=nky-1; m++){
00391 for (l=0; l<=nkz-1; l++){
00392 phi = 2.0 * pi * ran2(&rseed);
00393 phi2 = 2.0 * pi * ran2(&rseed);
00394 phi3 = 2.0 * pi * ran2(&rseed);
00395 amp = (1.0/sqrt(kx[n]*kx[n] + ky[m]*ky[m] + kz[l]*kz[l]));
00396 for (k=ks; k<=ke+1; k++) {
00397 za=(float)(k-ks)*(1.0/(float)(ke-ks+1));
00398 for (j=js; j<=je+1; j++) {
00399 ya=(float)(j-js)*(1.0/(float)(je-js+1));
00400 for (i=is; i<=ie+1; i++) {
00401 xa=(float)(i-is)*(2.0/(float)(ie-is+1));
00402 tmp = kx[n]*xa + ky[m]*ya + kz[l]*za;
00403 az[k][j][i] += amp*sin(tmp + phi);
00404 ay[k][j][i] += amp*sin(tmp + phi2);
00405 ax[k][j][i] += amp*sin(tmp + phi3);
00406 }
00407 }
00408 }
00409 }
00410 }
00411 }
00412
00413 printf("za=%f\n",za);
00414 Bsum = 0.0;
00415 #ifdef MHD
00416 dp = par_getd("problem","dp");
00417 for (k=ks; k<=ke; k++) {
00418 for (j=js; j<=je; j++) {
00419 for (i=is; i<=ie; i++) {
00420 pGrid->U[k][j][i].d = rho0;
00421 pGrid->U[k][j][i].M1 = 0.0;
00422 pGrid->U[k][j][i].M2 = 0.0;
00423 pGrid->U[k][j][i].M3 = 0.0;
00424 pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00425 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00426 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00427 pGrid->B1i[k][j][i] = (az[k][j+1][i] - az[k][j][i]) -
00428 (ay[k+1][j][i] - ay[k][j][i]);
00429 pGrid->B2i[k][j][i] = (ax[k+1][j][i] - ax[k][j][i]) -
00430 (az[k][j][i+1] - az[k][j][i]);
00431 pGrid->B3i[k][j][i] = (ay[k][j][i+1] - ay[k][j][i]) -
00432 (ax[k][j+1][i] - ax[k][j][i]);
00433
00434 Bsum += sqrt(pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00435 pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00436 }
00437 }
00438 }
00439 Btot = Bsum/((float)(je-js+1) * (float)(ie-is+1) * (float)(ke-ks+1));
00440 printf("Btot=%f\n",Btot);
00441 Bsum = 0.0;
00442 for (k=ks; k<=ke; k++) {
00443 for (j=js; j<=je; j++) {
00444 for (i=is; i<=ie; i++) {
00445 pGrid->B1i[k][j][i] = pGrid->B1i[k][j][i] * b0 / Btot;
00446 if (i == ie) pGrid->B1i[k][j][i+1] = pGrid->B1i[k][j][is];
00447 pGrid->B2i[k][j][i] = pGrid->B2i[k][j][i] * b0 / Btot;
00448 if (j == je) pGrid->B2i[k][j+1][i] = pGrid->B2i[k][js][i];
00449 pGrid->B3i[k][j][i] = pGrid->B3i[k][j][i] * b0 / Btot;
00450 if (k == ke) pGrid->B3i[k+1][j][i] = pGrid->B3i[ks][j][i];
00451
00452 Bsum += sqrt(pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00453 pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00454 }
00455 }
00456 }
00457
00458 for (k=ks; k<=ke; k++) {
00459 for (j=js; j<=je; j++) {
00460 for (i=is; i<=ie; i++) {
00461 pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00462 pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00463 pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00464 pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00465 + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00466 }
00467 }
00468 }
00469 Btot = Bsum/((float)(je-js+1) * (float)(ie-is+1) * (float)(ke-ks+1));
00470 printf("Btot/b0=%e %e\n",Btot,b0);
00471 #endif //MHD
00472 }
00473
00474 CoolingFunc = KoyInut;
00475 }
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00488 {
00489 return;
00490 }
00491
00492 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00493 {
00494 return;
00495 }
00496
00497
00498
00499 static Real logd(const Grid *pG, const int i, const int j, const int k)
00500 {
00501 return log10(pG->U[k][j][i].d);
00502 }
00503
00504 Gasfun_t get_usr_expr(const char *expr)
00505 {
00506 if(strcmp(expr,"logd")==0) return logd;
00507 return NULL;
00508 }
00509
00510 VGFunout_t get_usr_out_fun(const char *name){
00511 return NULL;
00512 }
00513
00514
00515 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00516 {
00517 return;
00518 }
00519
00520 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00521 {
00522 return;
00523 }
00524
00525
00526
00527
00528
00529 #define IM1 2147483563
00530 #define IM2 2147483399
00531 #define AM (1.0/IM1)
00532 #define IMM1 (IM1-1)
00533 #define IA1 40014
00534 #define IA2 40692
00535 #define IQ1 53668
00536 #define IQ2 52774
00537 #define IR1 12211
00538 #define IR2 3791
00539 #define NTAB 32
00540 #define NDIV (1+IMM1/NTAB)
00541 #define RNMX (1.0-DBL_EPSILON)
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556 double ran2(long int *idum){
00557 int j;
00558 long int k;
00559 static long int idum2=123456789;
00560 static long int iy=0;
00561 static long int iv[NTAB];
00562 double temp;
00563
00564 if (*idum <= 0) {
00565 if (-(*idum) < 1) *idum=1;
00566 else *idum = -(*idum);
00567 idum2=(*idum);
00568 for (j=NTAB+7;j>=0;j--) {
00569 k=(*idum)/IQ1;
00570 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00571 if (*idum < 0) *idum += IM1;
00572 if (j < NTAB) iv[j] = *idum;
00573 }
00574 iy=iv[0];
00575 }
00576 k=(*idum)/IQ1;
00577 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00578 if (*idum < 0) *idum += IM1;
00579 k=idum2/IQ2;
00580 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00581 if (idum2 < 0) idum2 += IM2;
00582 j=(int)(iy/NDIV);
00583 iy=iv[j]-idum2;
00584 iv[j] = *idum;
00585 if (iy < 1) iy += IMM1;
00586 if ((temp=AM*iy) > RNMX) return RNMX;
00587 else return temp;
00588 }
00589
00590 #undef IM1
00591 #undef IM2
00592 #undef AM
00593 #undef IMM1
00594 #undef IA1
00595 #undef IA2
00596 #undef IQ1
00597 #undef IQ2
00598 #undef IR1
00599 #undef IR2
00600 #undef NTAB
00601 #undef NDIV
00602 #undef RNMX