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