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 <string.h>
00030 #include <time.h>
00031 #include "defs.h"
00032 #include "athena.h"
00033 #include "globals.h"
00034 #include "prototypes.h"
00035 #include "particles/particle.h"
00036
00037 #ifndef SHEARING_BOX
00038 #error : The streaming3d problem requires shearing-box to be enabled.
00039 #endif
00040
00041 #ifndef PARTICLES
00042 #error : The streaming3d problem requires particles to be enabled.
00043 #endif
00044
00045 #ifndef ISOTHERMAL
00046 #error : The streaming3d problem requires isothermal equation of state.
00047 #endif
00048
00049
00050
00051 Real vsc1,vsc2;
00052 int ipert;
00053
00054 Real x1min,x1max,x2min,x2max,x3min,x3max,Lx,Ly,Lz,Lg;
00055 long Npar;
00056
00057 long ntrack;
00058 long nlis;
00059 int mytype;
00060 Real dpar_thresh;
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 double ran2(long int *idum);
00074 double Normal(long int *idum);
00075 Real Erf(Real z);
00076
00077 void MultiNSH(int n, Real *tstop, Real *mratio, Real etavk,
00078 Real *uxNSH, Real *uyNSH, Real *wxNSH, Real *wyNSH);
00079 static Real hst_rho_Vx_dVy(const Grid *pG,const int i,const int j,const int k);
00080 static void close_ix3(Grid *pGrid);
00081 static void close_ox3(Grid *pGrid);
00082 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3);
00083
00084 static int property_limit(const Grain *gr, const GrainAux *grsub);
00085 static int property_trace(const Grain *gr, const GrainAux *grsub);
00086 static int property_type(const Grain *gr, const GrainAux *grsub);
00087 static int property_dense(const Grain *gr, const GrainAux *grsub);
00088
00089 extern Real expr_dpar(const Grid *pG, const int i, const int j, const int k);
00090 extern Real expr_V1par(const Grid *pG, const int i, const int j, const int k);
00091 extern Real expr_V2par(const Grid *pG, const int i, const int j, const int k);
00092 extern Real expr_V3par(const Grid *pG, const int i, const int j, const int k);
00093 extern Real expr_V2(const Grid *pG, const int i, const int j, const int k);
00094
00095
00096
00097
00098
00099
00100 void problem(Grid *pGrid, Domain *pDomain)
00101 {
00102 int i,j,k,ks,pt,tsmode,vertical_bc;
00103 long p,q;
00104 Real ScaleHg,tsmin,tsmax,tscrit,amin,amax,Hparmin,Hparmax;
00105 Real *ep,*ScaleHpar,epsum,mratio,pwind,rhoaconv,etavk;
00106 Real *epsilon,*uxNSH,*uyNSH,**wxNSH,**wyNSH;
00107 Real rhog,h,x1,x2,x3,t,x1p,x2p,x3p,zmin,zmax,dx3_1,b;
00108 long int iseed = -1;
00109
00110 if (par_geti("grid","Nx3") == 1) {
00111 ath_error("[par_strat3d]: par_strat3d only works for 3D problem.\n");
00112 }
00113
00114 #ifdef MPI_PARALLEL
00115 if (par_geti("parallel","NGrid_x3") > 2) {
00116 ath_error("[par_strat2d]: Domain decomposition in x3 (z) is must be no more\
00117 than 2!\n");
00118 }
00119 #endif
00120
00121
00122 x1min = pGrid->x1_0 + (pGrid->is + pGrid->idisp)*pGrid->dx1;
00123 x1max = pGrid->x1_0 + (pGrid->ie + pGrid->idisp + 1.0)*pGrid->dx1;
00124 Lx = x1max - x1min;
00125
00126 x2min = pGrid->x2_0 + (pGrid->js + pGrid->jdisp)*pGrid->dx2;
00127 x2max = pGrid->x2_0 + (pGrid->je + pGrid->jdisp + 1.0)*pGrid->dx2;
00128 Ly = x2max - x2min;
00129
00130 x3min = par_getd("grid","x3min");
00131 x3max = par_getd("grid","x3max");
00132 Lz = x3max - x3min;
00133
00134 Lg = nghost*pGrid->dx3;
00135
00136 ks = pGrid->ks;
00137
00138
00139 Omega_0 = par_getd("problem","omega");
00140 qshear = par_getd_def("problem","qshear",1.5);
00141 ipert = par_geti_def("problem","ipert",1);
00142 vsc1 = par_getd_def("problem","vsc1",0.05);
00143 vsc2 = par_getd_def("problem","vsc2",0.0);
00144
00145 vsc1 = vsc1 * Iso_csound;
00146 vsc2 = vsc2 * Iso_csound;
00147
00148 ScaleHg = Iso_csound/Omega_0;
00149
00150
00151 Npar = (long)(par_geti("particle","parnumgrid"));
00152
00153 pGrid->nparticle = Npar*pGrid->partypes;
00154 for (i=0; i<pGrid->partypes; i++)
00155 pGrid->grproperty[i].num = Npar;
00156
00157 if (pGrid->nparticle+2 > pGrid->arrsize)
00158 particle_realloc(pGrid, pGrid->nparticle+2);
00159
00160 ep = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00161 ScaleHpar = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00162
00163 epsilon = (Real*)calloc_1d_array(pGrid->partypes, sizeof(Real));
00164 wxNSH = (Real**)calloc_2d_array(pGrid->Nx3+1, pGrid->partypes,sizeof(Real));
00165 wyNSH = (Real**)calloc_2d_array(pGrid->Nx3+1, pGrid->partypes,sizeof(Real));
00166 uxNSH = (Real*)calloc_1d_array(pGrid->Nx3+1, sizeof(Real));
00167 uyNSH = (Real*)calloc_1d_array(pGrid->Nx3+1, sizeof(Real));
00168
00169
00170 tsmode = par_geti("particle","tsmode");
00171 if (tsmode == 3) {
00172 tsmin = par_getd("problem","tsmin");
00173 tsmax = par_getd("problem","tsmax");
00174 tscrit= par_getd("problem","tscrit");
00175
00176 for (i=0; i<pGrid->partypes; i++) {
00177 tstop0[i] = tsmin*exp(i*log(tsmax/tsmin)/MAX(pGrid->partypes-1,1.0));
00178 pGrid->grproperty[i].rad = tstop0[i];
00179
00180 if (tstop0[i] < tscrit) pGrid->grproperty[i].integrator = 3;
00181 }
00182 }
00183 else {
00184 amin = par_getd("problem","amin");
00185 amax = par_getd("problem","amax");
00186
00187 for (i=0; i<pGrid->partypes; i++)
00188 pGrid->grproperty[i].rad = amin*
00189 exp(i*log(amax/amin)/MAX(pGrid->partypes-1,1.0));
00190
00191 if (tsmode >= 2) {
00192
00193 rhoaconv = par_getd_def("problem","rhoaconv",1.0);
00194
00195 for (i=0; i<pGrid->partypes; i++)
00196 grrhoa[i]=pGrid->grproperty[i].rad*rhoaconv;
00197 }
00198
00199 if (tsmode == 3)
00200 alamcoeff = par_getd("problem","alamcoeff");
00201 }
00202
00203
00204 Hparmax = par_getd("problem","hparmax");
00205 Hparmin = par_getd("problem","hparmin");
00206 for (i=0; i<pGrid->partypes; i++)
00207 ScaleHpar[i] = Hparmax*
00208 exp(-i*log(Hparmax/Hparmin)/MAX(pGrid->partypes-1,1.0));
00209
00210 #ifdef FEEDBACK
00211 mratio = par_getd_def("problem","mratio",0.0);
00212 pwind = par_getd_def("problem","pwind",0.0);
00213 if (mratio < 0.0)
00214 ath_error("[par_strat2d]: mratio must be positive!\n");
00215
00216 epsum = 0.0;
00217 for (i=0; i<pGrid->partypes; i++)
00218 {
00219 ep[i] = pow(pGrid->grproperty[i].rad,pwind); epsum += ep[i];
00220 }
00221
00222 for (i=0; i<pGrid->partypes; i++)
00223 {
00224 ep[i] = mratio*ep[i]/epsum;
00225 pGrid->grproperty[i].m = sqrt(2.0*PI)*ScaleHg/Lz*ep[i]*
00226 pGrid->Nx1*pGrid->Nx2*pGrid->Nx3/Npar;
00227 }
00228 #else
00229 mratio = 0.0;
00230 for (i=0; i<pGrid->partypes; i++)
00231 ep[i] = 0.0;
00232 #endif
00233
00234
00235 for (k=pGrid->ks; k<=pGrid->ke+1; k++) {
00236
00237 h = pGrid->x3_0 + (k+pGrid->kdisp)*pGrid->dx3;
00238 q = k - ks;
00239 etavk = fabs(vsc1+vsc2*SQR(h));
00240
00241 for (i=0; i<pGrid->partypes; i++) {
00242 epsilon[i] = ep[i]/ScaleHpar[i]*exp(-0.5*SQR(h/ScaleHg)*(SQR(1.0/ScaleHpar[i])-1.0))/erf(Lz/(sqrt(8.0)*ScaleHpar[i]*ScaleHg));
00243
00244 if (tsmode != 3)
00245 tstop0[i] = get_ts(pGrid,i,exp(-0.5*SQR(h/ScaleHg)),Iso_csound,etavk);
00246 }
00247
00248 MultiNSH(pGrid->partypes, tstop0, epsilon, etavk,
00249 &uxNSH[q], &uyNSH[q], wxNSH[q], wyNSH[q]);
00250 }
00251
00252
00253 for (k=pGrid->ks; k<=pGrid->ke; k++) {
00254 for (j=pGrid->js; j<=pGrid->je; j++) {
00255 for (i=pGrid->is; i<=pGrid->ie; i++) {
00256 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00257
00258 rhog = exp(-0.5*SQR(x3/ScaleHg));
00259 pGrid->U[k][j][i].d = rhog;
00260
00261 if (ipert != 1) {
00262 pGrid->U[k][j][i].M1 = 0.5*rhog*(uxNSH[k-ks]+uxNSH[k-ks+1]);
00263 pGrid->U[k][j][i].M2 = 0.5*rhog*(uyNSH[k-ks]+uyNSH[k-ks+1]);
00264 } else {
00265 pGrid->U[k][j][i].M1 = 0.0;
00266 pGrid->U[k][j][i].M2 = 0.0;
00267 }
00268
00269 pGrid->U[k][j][i].M3 = 0.0;
00270 #ifndef FARGO
00271 pGrid->U[k][j][i].M2 -= qshear*rhog*Omega*x1;
00272 #endif
00273
00274 }}}
00275
00276
00277 p = 0;
00278 dx3_1 = 1.0/pGrid->dx3;
00279 zmin = pGrid->x3_0 + (pGrid->ks+pGrid->kdisp)*pGrid->dx3;
00280 zmax = pGrid->x3_0 + (pGrid->ke+pGrid->kdisp+1)*pGrid->dx3;
00281
00282 for (q=0; q<Npar; q++) {
00283
00284 for (pt=0; pt<pGrid->partypes; pt++) {
00285
00286 x1p = x1min + Lx*ran2(&iseed);
00287 x2p = x2min + Ly*ran2(&iseed);
00288 x3p = ScaleHpar[pt]*ScaleHg*Normal(&iseed);
00289 while ((x3p >= zmax) || (x3p < zmin))
00290 x3p = ScaleHpar[pt]*ScaleHg*Normal(&iseed);
00291
00292 pGrid->particle[p].property = pt;
00293 pGrid->particle[p].x1 = x1p;
00294 pGrid->particle[p].x2 = x2p;
00295 pGrid->particle[p].x3 = x3p;
00296
00297 if (ipert != 1) {
00298
00299 cellk(pGrid, x3p, dx3_1, &k, &b);
00300 k = k-pGrid->ks; b = b - pGrid->ks;
00301
00302 pGrid->particle[p].v1 = (k+1-b)*wxNSH[k][pt]+(b-k)*wxNSH[k+1][pt];
00303 pGrid->particle[p].v2 = (k+1-b)*wyNSH[k][pt]+(b-k)*wyNSH[k+1][pt];
00304
00305 } else {
00306
00307 pGrid->particle[p].v1 = 0.0;
00308 pGrid->particle[p].v2 = vsc1+vsc2*SQR(x2p);
00309
00310 }
00311
00312 pGrid->particle[p].v3 = 0.0;
00313 #ifndef FARGO
00314 pGrid->particle[p].v2 -= qshear*Omega*x1p;
00315 #endif
00316
00317 pGrid->particle[p].pos = 1;
00318 pGrid->particle[p].my_id = p;
00319 #ifdef MPI_PARALLEL
00320 pGrid->particle[p].init_id = pGrid->my_id;
00321 #endif
00322 p++;
00323 }}
00324
00325
00326 StaticGravPot = ShearingBoxPot;
00327
00328 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00329
00330
00331
00332
00333 nlis = par_geti_def("problem","nlis",pGrid->Nx1*pGrid->Nx2*pGrid->Nx3);
00334
00335
00336 ntrack = par_geti_def("problem","ntrack",2000);
00337
00338
00339 dpar_thresh = par_geti_def("problem","dpar_thresh",10.0);
00340
00341
00342 vertical_bc = par_geti_def("problem","vertical_bc",0);
00343
00344 if (vertical_bc == 1)
00345 {
00346 set_bvals_mhd_fun(left_x3, close_ix3);
00347 set_bvals_mhd_fun(right_x3, close_ox3);
00348 }
00349
00350
00351 free(ep); free(ScaleHpar);
00352 free(epsilon);
00353 free_2d_array(wxNSH); free_2d_array(wyNSH);
00354 free(uxNSH); free(uyNSH);
00355
00356 return;
00357 }
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00370 {
00371 return;
00372 }
00373
00374 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00375 {
00376
00377 StaticGravPot = ShearingBoxPot;
00378
00379 Omega_0 = par_getd("problem","omega");
00380 qshear = par_getd_def("problem","qshear",1.5);
00381 ipert = par_geti_def("problem","ipert",1);
00382
00383 x1min = pG->x1_0 + (pG->is + pG->idisp)*pG->dx1;
00384 x1max = pG->x1_0 + (pG->ie + pG->idisp + 1.0)*pG->dx1;
00385 Lx = x1max - x1min;
00386
00387 x2min = pG->x2_0 + (pG->js + pG->jdisp)*pG->dx2;
00388 x2max = pG->x2_0 + (pG->je + pG->jdisp + 1.0)*pG->dx2;
00389 Ly = x2max - x2min;
00390
00391 x3min = par_getd("grid","x3min");
00392 x3max = par_getd("grid","x3max");
00393 Lz = x3max - x3min;
00394
00395 Lg = nghost*pG->dx3;
00396
00397 vsc1 = par_getd_def("problem","vsc1",0.05);
00398 vsc2 = par_getd_def("problem","vsc2",0.0);
00399
00400 vsc1 = vsc1 * Iso_csound;
00401 vsc2 = vsc2 * Iso_csound;
00402
00403 Npar = (int)(sqrt(par_geti("particle","parnumgrid")));
00404 nlis = par_geti_def("problem","nlis",pG->Nx1*pG->Nx2*pG->Nx3);
00405 ntrack = par_geti_def("problem","ntrack",2000);
00406
00407 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00408
00409 return;
00410 }
00411
00412 static Real hst_rho_Vx_dVy(const Grid *pG,
00413 const int i, const int j, const int k)
00414 {
00415 Real x1,x2,x3;
00416 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00417 #ifdef FARGO
00418 return pG->U[k][j][i].M1*pG->U[k][j][i].M2/pG->U[k][j][i].d;
00419 #else
00420 return pG->U[k][j][i].M1*(pG->U[k][j][i].M2/pG->U[k][j][i].d
00421 + qshear*Omega_0*x1);
00422 #endif
00423 }
00424
00425 static void close_ix3(Grid *pGrid)
00426 {
00427 int ks = pGrid->ks;
00428 int i,j,k,il,iu,jl,ju;
00429
00430 iu = pGrid->ie + nghost;
00431 il = pGrid->is - nghost;
00432
00433 ju = pGrid->je + nghost;
00434 jl = pGrid->js - nghost;
00435
00436 for (k=1; k<=nghost; k++)
00437 for (j=jl; j<=ju; j++)
00438 for (i=il; i<=iu; i++) {
00439 pGrid->U[ks-k][j][i] = pGrid->U[ks][j][i];
00440 pGrid->U[ks-k][j][i].M3 = 0.0;
00441 }
00442
00443 return;
00444 }
00445
00446 static void close_ox3(Grid *pGrid)
00447 {
00448 int ke = pGrid->ke;
00449 int i,j,k,il,iu,jl,ju;
00450
00451 iu = pGrid->ie + nghost;
00452 il = pGrid->is - nghost;
00453
00454 ju = pGrid->je + nghost;
00455 jl = pGrid->js - nghost;
00456
00457 for (k=1; k<=nghost; k++)
00458 for (j=jl; j<=ju; j++)
00459 for (i=il; i<=iu; i++) {
00460 pGrid->U[ke+k][j][i] = pGrid->U[ke][j][i];
00461 pGrid->U[ke+k][j][i].M3 = 0.0;
00462 }
00463
00464 return;
00465 }
00466
00467 Gasfun_t get_usr_expr(const char *expr)
00468 {
00469 return NULL;
00470 }
00471
00472 VGFunout_t get_usr_out_fun(const char *name){
00473 return NULL;
00474 }
00475
00476 #ifdef PARTICLES
00477 PropFun_t get_usr_par_prop(const char *name)
00478 {
00479 if (strcmp(name,"limit")==0) return property_limit;
00480 if (strcmp(name,"trace")==0) return property_trace;
00481 if (strcmp(name,"type")==0) return property_type;
00482 if (strcmp(name,"dense")==0) return property_dense;
00483
00484 return NULL;
00485 }
00486
00487 void gasvshift(const Real x1, const Real x2, const Real x3,
00488 Real *u1, Real *u2, Real *u3)
00489 {
00490 return;
00491 }
00492
00493 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00494 const Real v1, const Real v2, const Real v3)
00495 {
00496 Real z,fac;
00497
00498 ft->x1 -= 2.0*(vsc1 + vsc2*SQR(x3))*Omega_0;
00499
00500 if(x3 > x3max)
00501 z = x3-Lz;
00502 else if (x3 < x3min)
00503 z = x3+Lz;
00504 else
00505 z = x3;
00506
00507 fac = Lg/(0.5*Lz+Lg-fabs(z));
00508 ft->x3 -= SQR(Omega_0)*z*(1.0-SQR(fac)*fac);
00509
00510
00511 return;
00512 }
00513 #endif
00514
00515 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00516 {
00517 return;
00518 }
00519
00520
00521
00522
00523
00524
00525
00526 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00527 {
00528 return;
00529 }
00530
00531
00532
00533
00534
00535
00536 static Real ShearingBoxPot(const Real x1, const Real x2, const Real x3)
00537 {
00538 Real z,z0,phi=0.0;
00539
00540 #ifndef FARGO
00541 phi -= qshear*SQR(Omega_0*x1);
00542 #endif
00543
00544 if(x3 > x3max)
00545 z=fabs(x3-Lz);
00546 else if (x3 < x3min)
00547 z=fabs(x3+Lz);
00548 else
00549 z=fabs(x3);
00550
00551 phi += 0.5*SQR(Omega_0*z);
00552
00553
00554 z0 = 0.5*Lz+Lg;
00555 phi -= SQR(Omega_0*Lg)*Lg*(2.0*z-z0)/(2.0*SQR(z0-z));
00556
00557
00558
00559 return phi;
00560 }
00561
00562
00563
00564 static int property_limit(const Grain *gr, const GrainAux *grsub)
00565 {
00566 if ((gr->pos == 1) && (gr->my_id<nlis))
00567 return 1;
00568 else
00569 return 0;
00570 }
00571
00572
00573
00574 static int property_trace(const Grain *gr, const GrainAux *grsub)
00575 {
00576 if ((gr->pos == 1) && (gr->my_id<ntrack))
00577 return 1;
00578 else
00579 return 0;
00580 }
00581
00582
00583
00584 static int property_type(const Grain *gr, const GrainAux *grsub)
00585 {
00586 if ((gr->pos == 1) && (gr->property == mytype))
00587 return 1;
00588 else
00589 return 0;
00590 }
00591
00592
00593
00594 static int property_dense(const Grain *gr, const GrainAux *grsub)
00595 {
00596 if ((gr->pos == 1) && (grsub->dpar > dpar_thresh))
00597 return 1;
00598 else
00599 return 0;
00600 }
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612 void MultiNSH(int n, Real *tstop, Real *mratio, Real etavk,
00613 Real *uxNSH, Real *uyNSH, Real *wxNSH, Real *wyNSH)
00614 {
00615 int i,j;
00616 Real *Lambda1,**Lam1GamP1, **A, **B, **Tmp;
00617
00618 Lambda1 = (Real*)calloc_1d_array(n, sizeof(Real));
00619 Lam1GamP1=(Real**)calloc_2d_array(n, n, sizeof(Real));
00620 A = (Real**)calloc_2d_array(n, n, sizeof(Real));
00621 B = (Real**)calloc_2d_array(n, n, sizeof(Real));
00622 Tmp = (Real**)calloc_2d_array(n, n, sizeof(Real));
00623
00624
00625 for (i=0; i<n; i++){
00626 for (j=0; j<n; j++)
00627 Lam1GamP1[i][j] = mratio[j];
00628 Lam1GamP1[i][i] += 1.0;
00629 Lambda1[i] = 1.0/(tstop[i]+1.0e-16);
00630 for (j=0; j<n; j++)
00631 Lam1GamP1[i][j] *= Lambda1[i];
00632 }
00633
00634
00635 MatrixMult(Lam1GamP1, Lam1GamP1, n,n,n, Tmp);
00636 for (i=0; i<n; i++) Tmp[i][i] += 1.0;
00637 InverseMatrix(Tmp, n, B);
00638 for (i=0; i<n; i++)
00639 for (j=0; j<n; j++)
00640 B[i][j] *= Lambda1[j];
00641 MatrixMult(Lam1GamP1, B, n,n,n, A);
00642
00643
00644 *uxNSH = 0.0; *uyNSH = 0.0;
00645 for (i=0; i<n; i++){
00646 wxNSH[i] = 0.0;
00647 wyNSH[i] = 0.0;
00648 for (j=0; j<n; j++){
00649 wxNSH[i] -= B[i][j];
00650 wyNSH[i] -= A[i][j];
00651 }
00652 wxNSH[i] *= 2.0*etavk;
00653 wyNSH[i] *= etavk;
00654 *uxNSH -= mratio[i]*wxNSH[i];
00655 *uyNSH -= mratio[i]*wyNSH[i];
00656 wyNSH[i] += etavk;
00657 }
00658
00659 free(Lambda1);
00660 free_2d_array(A); free_2d_array(B);
00661 free_2d_array(Lam1GamP1); free_2d_array(Tmp);
00662
00663 return;
00664 }
00665
00666
00667
00668
00669 #define IM1 2147483563
00670 #define IM2 2147483399
00671 #define AM (1.0/IM1)
00672 #define IMM1 (IM1-1)
00673 #define IA1 40014
00674 #define IA2 40692
00675 #define IQ1 53668
00676 #define IQ2 52774
00677 #define IR1 12211
00678 #define IR2 3791
00679 #define NTAB 32
00680 #define NDIV (1+IMM1/NTAB)
00681 #define RNMX (1.0-DBL_EPSILON)
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 double ran2(long int *idum)
00697 {
00698 int j;
00699 long int k;
00700 static long int idum2=123456789;
00701 static long int iy=0;
00702 static long int iv[NTAB];
00703 double temp;
00704
00705 if (*idum <= 0) {
00706 if (-(*idum) < 1) *idum=1;
00707 else *idum = -(*idum);
00708 idum2=(*idum);
00709 for (j=NTAB+7;j>=0;j--) {
00710 k=(*idum)/IQ1;
00711 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00712 if (*idum < 0) *idum += IM1;
00713 if (j < NTAB) iv[j] = *idum;
00714 }
00715 iy=iv[0];
00716 }
00717 k=(*idum)/IQ1;
00718 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00719 if (*idum < 0) *idum += IM1;
00720 k=idum2/IQ2;
00721 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00722 if (idum2 < 0) idum2 += IM2;
00723 j=(int)(iy/NDIV);
00724 iy=iv[j]-idum2;
00725 iv[j] = *idum;
00726 if (iy < 1) iy += IMM1;
00727 if ((temp=AM*iy) > RNMX) return RNMX;
00728 else return temp;
00729 }
00730
00731
00732
00733
00734 double Normal(long int *idum)
00735 {
00736 double Y,X1,X2;
00737
00738 X1 = ran2(idum);
00739 X2 = ran2(idum);
00740
00741 Y = sqrt(-2.0*log(X1+TINY_NUMBER))*cos(2*PI*X2);
00742
00743 return Y;
00744 }
00745
00746
00747
00748
00749 Real Erf(Real z)
00750 {
00751
00752 static double x[101]={
00753 0.000000e+000, 3.783387e-003, 7.709914e-003, 1.178500e-002, 1.601425e-002, 2.040352e-002,
00754 2.495885e-002, 2.968653e-002, 3.459307e-002, 3.968525e-002, 4.497008e-002, 5.045486e-002,
00755 5.614715e-002, 6.205480e-002, 6.818596e-002, 7.454909e-002, 8.115295e-002, 8.800667e-002,
00756 9.511969e-002, 1.025018e-001, 1.101632e-001, 1.181145e-001, 1.263667e-001, 1.349310e-001,
00757 1.438193e-001, 1.530440e-001, 1.626176e-001, 1.725534e-001, 1.828652e-001, 1.935671e-001,
00758 2.046738e-001, 2.162008e-001, 2.281639e-001, 2.405796e-001, 2.534651e-001, 2.668380e-001,
00759 2.807169e-001, 2.951209e-001, 3.100699e-001, 3.255844e-001, 3.416859e-001, 3.583966e-001,
00760 3.757395e-001, 3.937386e-001, 4.124186e-001, 4.318054e-001, 4.519256e-001, 4.728071e-001,
00761 4.944786e-001, 5.169701e-001, 5.403124e-001, 5.645379e-001, 5.896800e-001, 6.157732e-001,
00762 6.428537e-001, 6.709587e-001, 7.001271e-001, 7.303990e-001, 7.618162e-001, 7.944220e-001,
00763 8.282614e-001, 8.633812e-001, 8.998296e-001, 9.376570e-001, 9.769156e-001, 1.017659e+000,
00764 1.059945e+000, 1.103830e+000, 1.149376e+000, 1.196644e+000, 1.245701e+000, 1.296614e+000,
00765 1.349454e+000, 1.404292e+000, 1.461205e+000, 1.520272e+000, 1.581573e+000, 1.645193e+000,
00766 1.711221e+000, 1.779746e+000, 1.850864e+000, 1.924673e+000, 2.001274e+000, 2.080774e+000,
00767 2.163281e+000, 2.248909e+000, 2.337778e+000, 2.430008e+000, 2.525728e+000, 2.625070e+000,
00768 2.728170e+000, 2.835170e+000, 2.946219e+000, 3.061469e+000, 3.181080e+000, 3.305216e+000,
00769 3.434048e+000, 3.567755e+000, 3.706521e+000, 3.850536e+000, 4.000000e+000
00770 };
00771 static double y[101]={
00772 0.00000000e+000, 4.26907434e-003, 8.69953340e-003, 1.32973284e-002, 1.80686067e-002, 2.30197153e-002,
00773 2.81572033e-002, 3.34878242e-002, 3.90185379e-002, 4.47565113e-002, 5.07091186e-002, 5.68839404e-002,
00774 6.32887618e-002, 6.99315688e-002, 7.68205444e-002, 8.39640613e-002, 9.13706742e-002, 9.90491090e-002,
00775 1.07008250e-001, 1.15257124e-001, 1.23804883e-001, 1.32660778e-001, 1.41834139e-001, 1.51334337e-001,
00776 1.61170754e-001, 1.71352743e-001, 1.81889576e-001, 1.92790394e-001, 2.04064148e-001, 2.15719527e-001,
00777 2.27764884e-001, 2.40208149e-001, 2.53056730e-001, 2.66317410e-001, 2.79996226e-001, 2.94098338e-001,
00778 3.08627885e-001, 3.23587825e-001, 3.38979770e-001, 3.54803790e-001, 3.71058224e-001, 3.87739454e-001,
00779 4.04841688e-001, 4.22356710e-001, 4.40273635e-001, 4.58578645e-001, 4.77254725e-001, 4.96281391e-001,
00780 5.15634428e-001, 5.35285634e-001, 5.55202571e-001, 5.75348359e-001, 5.95681482e-001, 6.16155658e-001,
00781 6.36719759e-001, 6.57317799e-001, 6.77889021e-001, 6.98368078e-001, 7.18685336e-001, 7.38767318e-001,
00782 7.58537287e-001, 7.77916009e-001, 7.96822665e-001, 8.15175962e-001, 8.32895397e-001, 8.49902691e-001,
00783 8.66123358e-001, 8.81488386e-001, 8.95935967e-001, 9.09413237e-001, 9.21877939e-001, 9.33299942e-001,
00784 9.43662512e-001, 9.52963249e-001, 9.61214608e-001, 9.68443923e-001, 9.74692870e-001, 9.80016358e-001,
00785 9.84480847e-001, 9.88162149e-001, 9.91142807e-001, 9.93509200e-001, 9.95348535e-001, 9.96745927e-001,
00786 9.97781755e-001, 9.98529475e-001, 9.99054014e-001, 9.99410828e-001, 9.99645625e-001, 9.99794704e-001,
00787 9.99885782e-001, 9.99939162e-001, 9.99969080e-001, 9.99985060e-001, 9.99993164e-001, 9.99997050e-001,
00788 9.99998805e-001, 9.99999548e-001, 9.99999841e-001, 9.99999948e-001, 9.99999985e-001
00789 };
00790 double z1, g;
00791 int il, iu, im;
00792 z1 = fabs(z);
00793
00794 il=0; iu=100;
00795 while(iu-il>1) {
00796 im = (iu+il)/2;
00797 if (x[im] > z1) iu = im;
00798 else il = im;
00799 }
00800
00801 g = (y[iu]*(z1-x[il])+y[il]*(x[iu]-z1))/(x[iu]-x[il]);
00802
00803 g = MIN(g,1.0);
00804 if (z<0.0) g = -g;
00805 return g;
00806 }
00807
00808 #undef IM1
00809 #undef IM2
00810 #undef AM
00811 #undef IMM1
00812 #undef IA1
00813 #undef IA2
00814 #undef IQ1
00815 #undef IQ2
00816 #undef IR1
00817 #undef IR2
00818 #undef NTAB
00819 #undef NDIV
00820 #undef RNMX