00001
00002 c
00003 c
00004 c
00005 c
00006 c
00007 c
00008 c
00009 c
00010
00011 c Fortran-77 worker function (hidden inside mk_aniso_king.C).
00012
00013 subroutine aking(mass, xx, yy, zz, vx, vy, vz,
00014 & alpha1_in, alpha3_in, iseed, N_in, W0_in, debug,
00015 & potential, kinetic, tidal_energy,
00016 & coord)
00017
00018 c Program to make an anisotropic King model (cf. Heggie & Ramamani,
00019 c MNRAS, 272, 317-22, 1995.)
00020 c
00021 c Written (thrown together) by d.c.heggie@uk.ac.ed
00022 c This program comes without warranty.
00023 c
00024 c Function version created by SLWM.
00025 c
00026 c Arguments are: arrays for masses, coordinates and velocities,
00027 c alpha1 (< 0), alpha3, iseed, N, W0 (1 < W0 < 12),
00028 c debug.
00029 c
00030 c The Galactic tidal potential per unit mass (with the convention
00031 c that the potential is negative deep inside a bound system) is
00032 c
00033 c 0.5 * (alpha1*x**2 + alpha3*z**2)
00034 c
00035 c For a point-mass tidal field, we have
00036 c
00037 c alpha1 = -3 Omega^2
00038 c alpha3 = Omega^2
00039 c
00040 c where Omega^2 is as defined in Binney and Tremaine (and used in
00041 c kira).
00042 c
00043 c iseed is a seed for setting the random number generator
00044 c N is the number of particles required in an N-body model; if n<=1
00045 c no particles are generated, but the surface density profile is
00046 c computed and printed
00047 c W0 is the concentration of a King model from which the anisotropic
00048 c model is generated by perturbation
00049 c |debug| > 0 prints/saves info about the King model itself
00050 c debug < 0 prints extra diagnostic messages
00051 c
00052 c Parameters from original program now set to specific values:
00053 c
00054 c mtot is the total mass of the system
00055 c G is G
00056 c xstart generates the initial value of W by W(0) = W0 - xstart
00057 c
00058 c **** Take mtot = G = 1, xstart = 0.001. ****
00059 c
00060 c Units:
00061 c
00062 c "Internal" code length unit (used in computation of the King model)
00063 c is the core radius (defined in the usual way).
00064 c
00065 c "External" units are defined by the fact that G = 1, the model has
00066 c total mass mtot = 1, and exactly fills its Roche lobe in the Galactic
00067 c tidal field (i.e. defined by G, mtot, and the alpha parameters).
00068 c
00069 c The value of alpha1 determines the overall length scale of the output
00070 c model. If we plan to rescale externally (e.g. N-body units), then
00071 c this scale is irrelevant. ONLY -alpha3/alpha1, WHICH DEFINES THE
00072 c "SHAPE" OF THE TIDAL POTENTIAL, MATTERS IN THAT CASE.
00073 c
00074 c Output:
00075 c
00076 c Unit 0 (stderr) shows
00077 c 1. the course of the iteration in which the analytic
00078 c model is computed (debug < 0),
00079 c 2. the density profiles along the three axes (in
00080 c program units) (|debug| > 0), and,
00081 c 3. the limiting radii along the three axes, in both
00082 c program and user units
00083 c Unit 7 produces surface densities (along the x and z axes,
00084 c as viewed along the y-axis) if N<=1 (|debug| > 1)
00085 c Unit 8 produces a list of the analytic solution (|debug| > 1)
00086
00087 implicit double precision (a-h,o-z)
00088
00089 double precision mass(N_in), xx(N_in), yy(N_in), zz(N_in),
00090 & vx(N_in), vy(N_in), vz(N_in)
00091 double precision alpha1, alpha3, W0
00092 integer iseed, N_in, debug
00093 double precision kinetic, potential, tidal_energy
00094
00095 parameter(ne=5,m=51,nb=3,nci=ne,ncj=ne-nb+1,nck=m+1,nsi=ne,
00096 & nsj=2*ne+1,nyj=ne,nyk=m,nstarm=10000)
00097
00098 c Note: nstarm was used to set the size of the coord array
00099 c -- no longer used (dynamically allocated in mk_aniso_king).
00100
00101 common x(m),h,mm,n,c2,anorm
00102 common /path/ kmax,kount,dxsav,xp(200),yp(10,200)
00103 double precision mtot,mstar,j22,ke,kevir
00104 common /params/ w0,alpha1,alpha3,rho(m),rhodot(m),rlag,yking(2,m),
00105 & den,rking
00106 dimension scalv(ne),indexv(ne),y(ne,m),c(nci,ncj,nck),s(nsi,nsj)
00107 parameter (nmax=10)
00108
00109 dimension ystart(nmax),coord(N_in,3) ! <-- new declaration of coord
00110
00111 c data scalv/4*0.001d0,1.d0/
00112 data scalv/5*1.d0/
00113 external derivs,rkqc
00114
00115 c Function erf could/should be replaced by intrinsic if possible.
00116
00117 densty(z) = -sqrt(z)*(z+1.5d0)+0.75*sqrt(pi)*exp(z)
00118 & *dch_erf(sqrt(z))
00119
00120 if (alpha3_in.ne.0.d0 .and. alpha1_in.ge.0.d0) then
00121 write(0,*) 'aking: alpha1 should be negative'
00122 stop
00123 endif
00124 if (w0_in.lt.0.1d0) then
00125 write(0,*) 'aking: w0 too small'
00126 stop
00127 endif
00128 if (w0_in.gt.12.d0) then
00129 write(0,*) 'aking: w0 too large'
00130 stop
00131 endif
00132
00133 pi = 4.d0*atan(1.d0)
00134 itmax=100
00135 conv=5.e-6
00136 h=1.d0/(m-1)
00137
00138 mtot = 1.0
00139 G = 1.0
00140
00141 alpha1 = alpha1_in
00142 alpha3 = alpha3_in
00143 N = N_in
00144 W0 = W0_in
00145
00146 xstart = 0.001
00147
00148 itmax = 100
00149 c rmin = 0.2
00150 c xstart = 0.2
00151 slowc = 1.d0
00152 do i = 1,5
00153 indexv(i) = i
00154 enddo
00155 anorm=1.
00156 h = (w0-xstart)**2/(m-1)
00157 den = densty(w0)
00158 dxsav = 0.01d0*w0
00159 kmax = 200
00160 nvar = 2
00161 c xstart = 0.02
00162 c xstart = 0.001
00163 ystart(1) = 2.d0*xstart/3.d0
00164 ystart(2) = -2.d0/3.d0
00165
00166 c Compute the King model and perturbations to it.
00167 c
00168 c In the computation of the King model, y1 = r^2, y2 = d(y1)/dx,
00169 c and x is dimensionless depth (w).
00170
00171 x1 = w0 - xstart
00172 x2 = 0.d0
00173 tol = 1.d-6
00174 h1 = 0.01d0
00175 hmin = 0.0000000001d0
00176
00177 call odeint(ystart,nvar,x1,x2,tol,h1,hmin,nok,nbad,derivs,rkqc)
00178
00179 do k=1,m
00180 x(k) = w0-xstart-sqrt(h*(k-1))
00181 if (x(k).ge.0.d0) then
00182 rho(k) = densty(x(k))/den
00183 rhodot(k) = rho(k) + x(k)**1.5d0/den
00184 else
00185 rho(k) = 0.d0
00186 rhodot(k) = 0.d0
00187 endif
00188 if (x(k).gt.xp(1)) then
00189 yking(1,k) = (w0-x(k))*yp(1,1)/(w0 - xp(1))
00190 yking(2,k) = yp(2,1)
00191 else
00192 i = 1
00193 115 continue
00194 if ((x(k)-xp(i))*(x(k)-xp(i+1)).le.0.d0) then
00195 yking(1,k) = yp(1,i) + (yp(1,i+1)-yp(1,i))*
00196 & (x(k)-xp(i))/(xp(i+1)-xp(i))
00197 yking(2,k) = yp(2,i) + (yp(2,i+1)-yp(2,i))*
00198 & (x(k)-xp(i))/(xp(i+1)-xp(i))
00199 else
00200 i = i+1
00201 if (i.le.kount-1) then
00202 goto 115
00203 else
00204 if (debug .lt. 0) then
00205 write(0,*) 'aking: failing interpolation',
00206 & ' of King values?', xp(kount), x(k)
00207 endif
00208 yking(1,k) = yp(1,kount)
00209 yking(2,k) = yp(2,kount)
00210 endif
00211 endif
00212 endif
00213 y(1,k) = -0.000004*(1-x(k)/w0)
00214 y(2,k) = 0.000004/w0
00215 y(3,k) = 0.000001*(1-x(k)/w0)
00216 y(4,k) = -0.000001/w0
00217 y(5,k) = log(0.000002d0)
00218 enddo
00219
00220 call solvde(itmax,conv,slowc,scalv,indexv,ne,nb,m,y,nyj,nyk,
00221 & c,nci,ncj,nck,s,nsi,nsj,debug)
00222
00223 if (abs(debug) .gt. 1) then
00224 if (debug .lt. 0)
00225 $ write(0,*)'aking: writing model to unit 8...'
00226 do j = 1,m
00227 write(8,100) x(j),(yking(i,j),i=1,2),(y(i,j),i=1,5)
00228 100 format(0pf10.5,1p7e10.2)
00229 enddo
00230 endif
00231
00232 c Now compute the density profiles
00233
00234 s0 = 1.5d0*(1.d0+alpha3/alpha1) ! FIRST use of alpha3/alpha1
00235 s2x = -0.5d0*alpha3/alpha1 + 1.d0
00236 s2y = s2x - 1.5d0
00237 s2z = alpha3/alpha1 - 0.5d0
00238 eps = - exp(y(5,m))
00239 wyold = -1.d0
00240 wzold = -1.d0
00241
00242 if (abs(debug) .gt. 0)
00243 & write(0,*) ' radius density on x',
00244 & ' on y on z King'
00245
00246 fmax = 0.d0
00247 do i = 1,m
00248 r = sqrt(yking(1,i))
00249 wx = x(i) + y(1,i)*s0 + y(3,i)*s2x - eps*r**2*
00250 & (2.d0/3.d0)*(s0/3.d0 + s2x)
00251 wy = x(i) + y(1,i)*s0 + y(3,i)*s2y - eps*r**2*
00252 & (2.d0/3.d0)*(s0/3.d0 + s2y)
00253 wz = x(i) + y(1,i)*s0 + y(3,i)*s2z - eps*r**2*
00254 & (2.d0/3.d0)*(s0/3.d0 + s2z)
00255 if (wx.ge.0.d0) then
00256 rhox = densty(wx)/den
00257 else
00258 rhox = 0.d0
00259 endif
00260 if (wy.ge.0.d0) then
00261 rhoy = densty(wy)/den
00262 else
00263 rhoy = 0.d0
00264 endif
00265 if (wz.ge.0.d0) then
00266 rhoz = densty(wz)/den
00267 else
00268 rhoz = 0.d0
00269 endif
00270
00271 if (abs(debug) .gt. 0) write(0,101) r,rhox,rhoy,rhoz,rho(i)
00272 101 format(0pf10.5,1p4d15.5)
00273
00274 if (fmax.lt.r**2*rhox) fmax = r**2*rhox
00275 if (wyold.gt.0.d0) ymax = (r*wyold - wy*
00276 & sqrt(yking(1,i-1)))/(wyold-wy)
00277 if (wzold.gt.0.d0) zmax = (r*wzold - wz*
00278 & sqrt(yking(1,i-1)))/(wzold-wz)
00279 wyold = wy
00280 wzold = wz
00281 enddo
00282 r = r + 0.1
00283 10 continue
00284
00285 wc1 = -2.d0*yking(1,m)**1.5d0/yking(2,m)
00286
00287 c = 9 M(r) / (4 pi rho0)
00288
00289 wc2 = -wc1/sqrt(yking(1,m))
00290
00291 c = -9 M(r) / (4 pi rho0 r)
00292
00293 w0c1 = wc1*y(2,m)
00294 w0c2 = y(1,m) - w0c1/sqrt(yking(1,m))
00295 w2c1 = yking(1,m)**1.5d0*y(3,m)
00296 wx = wc1/r + wc2
00297 & +s0*(w0c1/r + w0c2) + s2x*w2c1/r**3 - eps*r**2*
00298 & (2.d0/3.d0)*(s0/3.d0 + s2x)
00299 wy = wc1/r + wc2
00300 & +s0*(w0c1/r + w0c2) + s2y*w2c1/r**3 - eps*r**2*
00301 & (2.d0/3.d0)*(s0/3.d0 + s2y)
00302 wz = wc1/r + wc2
00303 & +s0*(w0c1/r + w0c2) + s2z*w2c1/r**3 - eps*r**2*
00304 & (2.d0/3.d0)*(s0/3.d0 + s2z)
00305 if (wx.ge.0.d0) then
00306 rhox = densty(wx)/den
00307 else
00308 rhox = 0.d0
00309 endif
00310 if (wy.ge.0.d0) then
00311 rhoy = densty(wy)/den
00312 else
00313 rhoy = 0.d0
00314 endif
00315 if (wz.ge.0.d0) then
00316 rhoz = densty(wz)/den
00317 else
00318 rhoz = 0.d0
00319 endif
00320 if (wyold.gt.0.d0) ymax = (r*wyold - wy*(r-0.1d0))/
00321 & (wyold-wy)
00322 if (wzold.gt.0.d0) zmax = (r*wzold - wz*(r-0.1d0))/
00323 & (wzold-wz)
00324 wyold = wy
00325 wzold = wz
00326 rhok = 0.d0
00327
00328 if (abs(debug) .gt. 0) write(0,101) r,rhox,rhoy,rhoz,rhok
00329
00330 if (fmax.lt.r**2*rhox) fmax = r**2*rhox
00331 r = r + 0.1
00332 if (r.lt.rlag) goto 10
00333 fmax = 1.1d0
00334
00335 c Reminder from Steve to Steve (7/16/98): at this point, lengths
00336 c are measured in core radii and the density is scaled to the
00337 c central value.
00338
00339 write(0,*) 'aking: scaled masses: King ', wc1,
00340 & ' anisotropic ', wc1+s0*w0c1
00341 write(0,*) 'aking: radii (core radii):',rlag, ymax, zmax
00342 rking = sqrt(yking(1,m))
00343 write(0,*) 'aking: King cutoff radius (core radii):',rking
00344 write(0,'('' aking: scaled densities: King: '',1pe9.3, '
00345 &
00346 & wc1*3/(4*pi*rking**3),
00347 & (wc1 + s0*w0c1)*3/(4*pi*rlag**3), W0
00348
00349 c Next section computes surface density viewed parallel to y-axis
00350
00351 if (abs(debug).gt.1) then
00352 zero = 0.d0
00353 sigma0 = sigma(zero,zero,y)
00354 if (debug .lt. 0)
00355 & write(0,*) 'Writing surface densities to unit 7...'
00356 do i = 1,250
00357 x1 = 10.d0**(0.01d0*(i-1))
00358 if (x1.le.rlag) then
00359 sigmax = sigma(x1,zero,y)/sigma0
00360 sigmaz = sigma(zero,x1,y)/sigma0
00361 write(7,*) x1, sigmax, sigmaz
00362 endif
00363 enddo
00364 endif
00365
00366 c Next section computes N-body initial conditions. Note that alpha1
00367 c is used to establish scaling.
00368
00369 idum = iseed
00370
00371 c Compute core radius in user units:
00372
00373 rc3 = 2.d0*g*mtot*eps/(alpha1*(wc1 + s0*w0c1)) ! FIRST use of alpha1
00374
00375 c (This eps really is the eps in the paper, although it is
00376 c obtained here as a solution to a boundary-value problem.)
00377
00378 rc = rc3**(1.d0/3.d0)
00379 write(0,*)'aking: core radius (user units) = ', rc
00380
00381 c Compute 1 / (1-D velocity dispersion) in user units:
00382
00383 rc2 = rc**2
00384 j22 = 2.d0*eps/(rc2*alpha1)
00385
00386 pe = 0.d0
00387 ke = 0.d0
00388 sumx = 0.d0
00389 sumy = 0.d0
00390 sumz = 0.d0
00391
00392 do 40 i = 1,n
00393 if (n.le.1) goto 40
00394 15 continue
00395 rstar = rlag*ran2(idum)
00396 costh = 2.d0*ran2(idum)-1.d0
00397 phi = 2.d0*pi*ran2(idum)
00398 sinth = sqrt(1.d0-costh**2)
00399 xstar = rstar*sinth*cos(phi)
00400 ystar = rstar*sinth*sin(phi)
00401 zstar = rstar*costh
00402 r2s2 = (alpha3/alpha1-0.5d0)*(zstar**2-0.5d0*(xstar**2+
00403 & ystar**2))+ 0.75d0*(xstar**2-ystar**2)
00404 r2 = xstar**2 + ystar**2 + zstar**2
00405 s2 = r2s2/r2
00406 r = sqrt(r2)
00407 if (r.lt.rking) then
00408 if (r.lt.sqrt(yking(1,1))) then
00409 w1 = x(1) + y(1,1)*s0 + y(3,1)*s2 - eps*r**2*
00410 & (2.d0/3.d0)*(s0/3.d0 + s2)
00411 w = w0 - (r2/yking(1,1))*(w0 - w1)
00412 else
00413 j = 1
00414 20 continue
00415 j = j + 1
00416 if (r.gt.sqrt(yking(1,j))) goto 20
00417 wj1 = x(j-1) + y(1,j-1)*s0 + y(3,j-1)*s2 - eps*r**2*
00418 & (2.d0/3.d0)*(s0/3.d0 + s2)
00419 wj = x(j) + y(1,j)*s0 + y(3,j)*s2 - eps*r**2*
00420 & (2.d0/3.d0)*(s0/3.d0 + s2)
00421 w = wj1 + (r2-yking(1,j-1))*(wj-wj1)/
00422 & (yking(1,j)-yking(1,j-1))
00423 endif
00424 else
00425 w = wc1/r + wc2
00426 & +s0*(w0c1/r + w0c2) + s2*w2c1/r**3 - eps*r**2*
00427 & (2.d0/3.d0)*(s0/3.d0 + s2)
00428 endif
00429 if (w.gt.0.d0) then
00430 rhost = densty(w)/den
00431 else
00432 rhost = 0
00433 endif
00434 if (rstar**2*rhost.lt.fmax*ran2(idum)) go to 15
00435
00436 c Now choose speed
00437
00438 vmax = sqrt(2.d0*w)
00439 30 continue
00440 speed = vmax*ran2(idum)
00441 fstar = speed**2*(exp(-0.5d0*speed**2)-exp(-w))
00442 if (fstar.lt.2.d0*ran2(idum)/exp(1.d0)) go to 30
00443 costh = 2.d0*ran2(idum)-1.d0
00444 phi = 2.d0*pi*ran2(idum)
00445 sinth = sqrt(1.d0-costh**2)
00446
00447 c Convert positions and velocities from code units to user units.
00448
00449 speed = speed/sqrt(j22)
00450
00451 xstar = xstar*rc
00452 ystar = ystar*rc
00453 zstar = zstar*rc
00454 ustar = speed*sinth*cos(phi)
00455 vstar = speed*sinth*sin(phi)
00456 wstar = speed*costh
00457 mstar = mtot/n
00458
00459 coord(i,1) = xstar
00460 coord(i,2) = ystar
00461 coord(i,3) = zstar
00462
00463 c Ultimately, might be better to compute energies, etc. in
00464 c mk_aniso_king, but keep code as is for now. Avoid significant
00465 c cost of computing energy twice by passing energies back to
00466 c mk_aniso_king via the argument list.
00467
00468 if (i.gt.1) then
00469 do j = 1,i-1
00470 r2 = (coord(j,1)-xstar)**2 + (coord(j,2)-ystar)**2 +
00471 & (coord(j,3)-zstar)**2
00472 pe = pe - 1.d0/sqrt(r2)
00473 enddo
00474 endif
00475 ke = ke + mstar*speed**2
00476 sumx = sumx + mstar*xstar**2
00477 sumy = sumy + mstar*ystar**2
00478 sumz = sumz + mstar*zstar**2
00479
00480 mass(i) = mstar
00481 xx(i) = xstar
00482 yy(i) = ystar
00483 zz(i) = zstar
00484 vx(i) = ustar
00485 vy(i) = vstar
00486 vz(i) = wstar
00487 40 continue
00488 pe = g*pe*mstar**2
00489 e_tidal = 0.5*alpha1*sumx + 0.5*alpha3*sumz
00490 ke = ke*0.5d0
00491 xmax = rc*rlag
00492 ymax = rc*ymax
00493 zmax = rc*zmax
00494
00495 write(0,*) 'aking: radii (user units):', xmax, ymax, zmax
00496 write(0,*) 'aking: "Jacobi radius" (from alpha1 ',
00497 $ 'with phi ~ -1/r; user units):',
00498 $ (-g*mtot/alpha1)**(0.3333333333)
00499
00500
00501 if (n.gt.1) then
00502 rvir = -g*mtot**2/(2.d0*pe)
00503 rvir_tidal = -g*mtot**2/(2.d0*(pe+e_tidal))
00504 write(0,*) 'aking: virial radius of N-body system ',
00505 $ '(user units):', rvir
00506 write(0,*) 'aking: ....same, but with tidal potential ',
00507 $ 'included:', rvir_tidal
00508
00509 endif
00510
00511 c Derived quantities (user units) returned to calling program:
00512
00513 potential = pe
00514 kinetic = ke
00515 tidal_energy = e_tidal
00516
00517 return
00518 end
00519
00520 c
00521
00522 subroutine difeq(k,k1,k2,jsf,is1,isf,indexv,ne,s,nsi,nsj,y,nyj,nyk
00523 & )
00524 implicit double precision (a-h,o-z)
00525 parameter(m=51)
00526 common x(m),h,mm,n,c2,anorm
00527 common /params/ w0,alpha1,alpha3,rho(m),rhodot(m),rlag,
00528 & yking(2,m),den,rking
00529 dimension y(nyj,nyk),s(nsi,nsj),indexv(nyj)
00530
00531 g2(i) = -2.25d0*yking(2,i)**2*(rhodot(i)*(y(1,i) + (2.d0/9.d0)*
00532 & exp(y(5,i))*yking(1,i)) - y(2,i)*rho(i))/yking(1,i)
00533 g4(i) = 0.25d0*yking(2,i)**2*(6.d0*y(3,i)/yking(1,i)-
00534 & 9.d0*rhodot(i)*(y(3,i) + (2.d0/3.d0)
00535 & *exp(y(5,i))*yking(1,i))
00536 & + 9.d0*y(4,i)*rho(i))/yking(1,i)
00537
00538 if(k.eq.k1) then
00539 do i = 3,5
00540 do j = 6,10
00541 s(i,j) = 0.d0
00542 enddo
00543 enddo
00544 eps = -exp(y(5,1))
00545 b0 = 0.1d0*eps*rhodot(1)
00546 s(3,5+1) = 1.d0
00547 s(3,5+5) = -b0*yking(1,1)**2
00548 s(3,jsf) = y(1,1) - b0*yking(1,1)**2
00549 s(4,5+2) = 1.d0
00550 s(4,5+5) = -2.d0*b0*yking(1,1)*yking(2,1)
00551 s(4,jsf) = y(2,1) - 2.d0*b0*yking(1,1)*yking(2,1)
00552 s(5,5+3) = yking(2,1)
00553 s(5,5+4) = -yking(1,1)
00554 s(5,jsf) = y(3,1)*yking(2,1) - yking(1,1)*y(4,1)
00555 else if(k.gt.k2) then
00556 do i = 1,2
00557 do j = 6,10
00558 s(i,j) = 0.d0
00559 enddo
00560 enddo
00561 s(1,5+3) = 1.5d0*yking(2,m)
00562 s(1,5+4) = yking(1,m)
00563 s(1,jsf) = yking(1,m)*y(4,m) + 1.5d0*y(3,m)*yking(2,m)
00564 c now find radius where acceleration vanishes
00565 r = sqrt(yking(1,m))
00566 s0 = 1.5d0*(1.d0+alpha3/alpha1)
00567 s2 = -0.5d0*alpha3/alpha1 + 1.d0
00568 iter = 0
00569 eps = -exp(y(5,m))
00570 wc1 = -2.d0*yking(1,m)**1.5d0/yking(2,m)
00571 wc2 = -wc1/sqrt(yking(1,m))
00572 w0c1 = wc1*y(2,m)
00573 w0c2 = y(1,m) - w0c1/sqrt(yking(1,m))
00574 w2c1 = yking(1,m)**1.5d0*y(3,m)
00575 10 continue
00576 iter = iter + 1
00577 if (iter.gt.20) then
00578 write(0,*) 'aking: too many iterations finding r'
00579 write(0,*) (y(i,m),i=1,5)
00580 stop
00581 endif
00582 f = -wc1/r**2 +s0*(-w0c1/r**2)
00583 & -s2*3.d0*w2c1/r**4 - 2.d0*eps*r
00584 fdash = 2.d0*wc1/r**3 + 2.d0*s0*w0c1
00585 & /r**3 +12.d0*s2*w2c1/r**5 -
00586 & 2.d0*eps
00587 rnew = r - f/fdash
00588 if (abs((r-rnew)/rnew).gt.1.d-4) then
00589 r = rnew
00590 go to 10
00591 endif
00592 r = rnew
00593 rlag = r
00594 temp = 1.d0 - sqrt(yking(1,m))/r
00595 s(2,5+1) = s0
00596 s(2,5+2) = s0*2.d0*(yking(1,m)/yking(2,m))*temp
00597 s(2,5+3) = s2*yking(1,m)**1.5d0/r**3
00598 s(2,5+5) = exp(y(5,m))*r**2
00599 s(2,jsf) = 2.d0*yking(1,m)*temp/yking(2,m) +
00600 & s0*(y(1,m) + 2.d0*yking(1,m)*y(2,m)*temp/yking(2,m))
00601 & + s2*yking(1,m)**1.5d0*y(3,m)/r**3
00602 & - eps*r**2
00603 else
00604 do i = 1,5
00605 do j = 1,10
00606 s(i,j)=0.d0
00607 enddo
00608 enddo
00609 halfh = 0.5d0*(x(k)-x(k-1))
00610 s(1,1) = -1.d0
00611 s(1,2) = -halfh
00612 s(1,5+1)=1.d0
00613 s(1,5+2)=-halfh
00614 s(1,jsf)=y(1,k)-y(1,k-1)-halfh*(y(2,k)+y(2,k-1))
00615 s(2,1) = -halfh*(-2.25d0*(yking(2,k-1)**2/yking(1,k-1))*
00616 & rhodot(k-1))
00617 s(2,2) = -1.d0 - halfh*(-2.25d0*(yking(2,k-1)**2/
00618 & yking(1,k-1))*
00619 & (-rho(k-1)))
00620 s(2,5) = -halfh*(-2.25d0*(yking(2,k-1)**2/yking(1,k-1))*
00621 & rhodot(k-1)*
00622 & (2.d0/9.d0)*exp(y(5,k-1))*yking(1,k-1))
00623 s(2,5+1) = -halfh*(-2.25d0*(yking(2,k)**2/
00624 & yking(1,k))*rhodot(k))
00625 s(2,5+2) = 1.d0 - halfh*(-2.25d0*(yking(2,k)**2/yking(1,k))*
00626 & (-rho(k)))
00627 s(2,5+5) = -halfh*(-2.25d0*(yking(2,k)**2/
00628 & yking(1,k))*rhodot(k)*
00629 & (2.d0/9.d0)*exp(y(5,k))*yking(1,k))
00630 s(2,jsf) = y(2,k)-y(2,k-1)-halfh*(g2(k-1)+g2(k))
00631 s(3,3) = -1.d0
00632 s(3,4) = -halfh
00633 s(3,5+3)=1.d0
00634 s(3,5+4)=-halfh
00635 s(3,jsf)=y(3,k)-y(3,k-1)-halfh*(y(4,k)+y(4,k-1))
00636 s(4,3) = -halfh*(0.25d0*(yking(2,k-1)**2/
00637 & yking(1,k-1))*(6.d0/
00638 & yking(1,k-1)
00639 & -9.d0*rhodot(k-1)))
00640 s(4,4) = -1.d0 - halfh*(0.25d0*(yking(2,k-1)**2/
00641 & yking(1,k-1))*
00642 & 9.d0*
00643 & rho(k-1))
00644 s(4,5) = -halfh*(0.25d0*(yking(2,k-1)**2/
00645 & yking(1,k-1))*(-9.d0)*
00646 & rhodot(k-1)*(2.d0/3.d0)*exp(y(5,k-1))*yking(1,k-1))
00647 s(4,5+3) = -halfh*(0.25d0*(yking(2,k)**2/
00648 & yking(1,k))*(6.d0/
00649 & yking(1,k)
00650 & -9.d0*rhodot(k)))
00651 s(4,5+4) = 1.d0 - halfh*(0.25d0*(yking(2,k)**2/
00652 & yking(1,k))*9.d0*
00653 & rho(k))
00654 s(4,5+5) = -halfh*(0.25d0*(yking(2,k)**2/
00655 & yking(1,k))*(-9.d0)*
00656 & rhodot(k)*(2.d0/3.d0)*exp(y(5,k))*yking(1,k))
00657 s(4,jsf) = y(4,k)-y(4,k-1)-halfh*(g4(k-1)+g4(k))
00658 s(5,5) = -1.d0
00659 s(5,5+5) = 1.d0
00660 s(5,jsf) = -y(5,k-1) + y(5,k)
00661 endif
00662
00663 return
00664 end
00665
00666 c
00667
00668 subroutine solvde(itmax,conv,slowc,scalv,indexv,ne,nb,m,
00669 & y,nyj,nyk,c,nci,ncj,nck,s,nsi,nsj,debug)
00670 implicit double precision (a-h,o-z)
00671 integer debug
00672 parameter (nmax=10)
00673 dimension y(nyj,nyk),c(nci,ncj,nck),s(nsi,nsj),scalv(nyj),indexv(n
00674 & yj)
00675 dimension ermax(nmax),kmax(nmax)
00676
00677 k1=1
00678 k2=m
00679 nvars=ne*m
00680 j1=1
00681 j2=nb
00682 j3=nb+1
00683 j4=ne
00684 j5=j4+j1
00685 j6=j4+j2
00686 j7=j4+j3
00687 j8=j4+j4
00688 j9=j8+j1
00689 ic1=1
00690 ic2=ne-nb
00691 ic3=ic2+1
00692 ic4=ne
00693 jc1=1
00694 jcf=ic3
00695 do 16 it=1,itmax
00696 k=k1
00697 call difeq(k,k1,k2,j9,ic3,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk)
00698 call pinvs(ic3,ic4,j5,j9,jc1,k1,c,nci,ncj,nck,s,nsi,nsj)
00699 do 11 k=k1+1,k2
00700 kp=k-1
00701 call difeq(k,k1,k2,j9,ic1,ic4,indexv,ne,s,nsi,nsj,y,nyj,
00702 & nyk)
00703 call red(ic1,ic4,j1,j2,j3,j4,j9,ic3,jc1,jcf,kp,
00704 & c,nci,ncj,nck,s,nsi,nsj)
00705 call pinvs(ic1,ic4,j3,j9,jc1,k,c,nci,ncj,nck,s,nsi,nsj)
00706 11 continue
00707 k=k2+1
00708 call difeq(k,k1,k2,j9,ic1,ic2,indexv,ne,s,nsi,nsj,y,nyj,nyk)
00709 call red(ic1,ic2,j5,j6,j7,j8,j9,ic3,jc1,jcf,k2,
00710 & c,nci,ncj,nck,s,nsi,nsj)
00711 call pinvs(ic1,ic2,j7,j9,jcf,k2+1,c,nci,ncj,nck,s,nsi,nsj)
00712 call bksub(ne,nb,jcf,k1,k2,c,nci,ncj,nck)
00713 err=0.
00714 do 13 j=1,ne
00715 jv=indexv(j)
00716 ermax(j)=0.
00717 errj=0.
00718 kmax(j)=0
00719 vmax=0.
00720 do 12 k=k1,k2
00721 vz=abs(c(j,1,k))
00722 if(vz.gt.vmax) then
00723 vmax=vz
00724 km=k
00725 endif
00726 errj=errj+vz
00727 12 continue
00728 err=err+errj/scalv(jv)
00729 ermax(j)=c(j,1,km)/scalv(jv)
00730 kmax(j)=km
00731 13 continue
00732 err=err/nvars
00733 fac=slowc/max(slowc,err)
00734 do 15 jv=1,ne
00735 j=indexv(jv)
00736 do 14 k=k1,k2
00737 y(j,k)=y(j,k)-fac*c(jv,1,k)
00738 14 continue
00739 15 continue
00740 if (debug .lt. 0)
00741 & write(0,100) it,err,fac,(kmax(j),ermax(j),j=1,ne)
00742 if (err.lt.conv) return
00743 16 continue
00744 write(0,*) 'aking: itmax exceeded'
00745 stop
00746 100 format(1x,i4,2f12.6,(/5x,i5,f12.6))
00747 end
00748
00749
00750 subroutine bksub(ne,nb,jf,k1,k2,c,nci,ncj,nck)
00751 implicit double precision (a-h,o-z)
00752 dimension c(nci,ncj,nck)
00753 nbf=ne-nb
00754 do 13 k=k2,k1,-1
00755 kp=k+1
00756 do 12 j=1,nbf
00757 xx=c(j,jf,kp)
00758 do 11 i=1,ne
00759 c(i,jf,k)=c(i,jf,k)-c(i,j,k)*xx
00760 11 continue
00761 12 continue
00762 13 continue
00763 do 16 k=k1,k2
00764 kp=k+1
00765 do 14 i=1,nb
00766 c(i,1,k)=c(i+nbf,jf,k)
00767 14 continue
00768 do 15 i=1,nbf
00769 c(i+nb,1,k)=c(i,jf,kp)
00770 15 continue
00771 16 continue
00772
00773 return
00774 end
00775
00776 c
00777
00778 subroutine pinvs(ie1,ie2,je1,jsf,jc1,k,c,nci,ncj,nck,s,nsi,nsj)
00779 implicit double precision (a-h,o-z)
00780 parameter (zero=0.,one=1.,nmax=10)
00781 dimension c(nci,ncj,nck),s(nsi,nsj),pscl(nmax),indxr(nmax)
00782
00783 je2=je1+ie2-ie1
00784 js1=je2+1
00785 do 12 i=ie1,ie2
00786 big=zero
00787 do 11 j=je1,je2
00788 if(abs(s(i,j)).gt.big) big=abs(s(i,j))
00789 11 continue
00790 if(big.eq.zero) then
00791 write(0,*) 'aking: singular matrix, row all 0'
00792 stop
00793 endif
00794 pscl(i)=one/big
00795 indxr(i)=0
00796 12 continue
00797 do 18 id=ie1,ie2
00798 piv=zero
00799 do 14 i=ie1,ie2
00800 if(indxr(i).eq.0) then
00801 big=zero
00802 do 13 j=je1,je2
00803 if(abs(s(i,j)).gt.big) then
00804 jp=j
00805 big=abs(s(i,j))
00806 endif
00807 13 continue
00808 if(big*pscl(i).gt.piv) then
00809 ipiv=i
00810 jpiv=jp
00811 piv=big*pscl(i)
00812 endif
00813 endif
00814 14 continue
00815 if(s(ipiv,jpiv).eq.zero) then
00816 write(0,*) 'aking: singular matrix'
00817 stop
00818 endif
00819 indxr(ipiv)=jpiv
00820 pivinv=one/s(ipiv,jpiv)
00821 do 15 j=je1,jsf
00822 s(ipiv,j)=s(ipiv,j)*pivinv
00823 15 continue
00824 s(ipiv,jpiv)=one
00825 do 17 i=ie1,ie2
00826 if(indxr(i).ne.jpiv) then
00827 if(s(i,jpiv).ne.zero) then
00828 dum=s(i,jpiv)
00829 do 16 j=je1,jsf
00830 s(i,j)=s(i,j)-dum*s(ipiv,j)
00831 16 continue
00832 s(i,jpiv)=zero
00833 endif
00834 endif
00835 17 continue
00836 18 continue
00837 jcoff=jc1-js1
00838 icoff=ie1-je1
00839 do 21 i=ie1,ie2
00840 irow=indxr(i)+icoff
00841 do 19 j=js1,jsf
00842 c(irow,j+jcoff,k)=s(i,j)
00843 19 continue
00844 21 continue
00845
00846 return
00847 end
00848
00849
00850 subroutine red(iz1,iz2,jz1,jz2,jm1,jm2,jmf,ic1,jc1,jcf,kc,
00851 & c,nci,ncj,nck,s,nsi,nsj)
00852 implicit double precision (a-h,o-z)
00853 dimension c(nci,ncj,nck),s(nsi,nsj)
00854
00855 loff=jc1-jm1
00856 ic=ic1
00857 do 14 j=jz1,jz2
00858 do 12 l=jm1,jm2
00859 vx=c(ic,l+loff,kc)
00860 do 11 i=iz1,iz2
00861 s(i,l)=s(i,l)-s(i,j)*vx
00862 11 continue
00863 12 continue
00864 vx=c(ic,jcf,kc)
00865 do 13 i=iz1,iz2
00866 s(i,jmf)=s(i,jmf)-s(i,j)*vx
00867 13 continue
00868 ic=ic+1
00869 14 continue
00870
00871 return
00872 end
00873
00874
00875 double precision function dch_erf(x)
00876 implicit double precision (a-h,o-z)
00877
00878 half = 0.5d0
00879 if(x.lt.0.)then
00880 dch_erf=-gammp(half,x**2)
00881 else
00882 dch_erf=gammp(half,x**2)
00883 endif
00884
00885 return
00886 end
00887
00888
00889 double precision function gammp(a,x)
00890 implicit double precision (a-h,o-z)
00891
00892 if(x.lt.0..or.a.le.0.) stop
00893 if(x.lt.a+1.)then
00894 call gser(gammp,a,x,gln)
00895 else
00896 call gcf(gammcf,a,x,gln)
00897 gammp=1.-gammcf
00898 endif
00899
00900 return
00901 end
00902
00903
00904 subroutine gser(gamser,a,x,gln)
00905 implicit double precision (a-h,o-z)
00906 parameter (itmax=100,eps=3.e-7)
00907
00908 gln=gammln(a)
00909 if(x.le.0.)then
00910 if(x.lt.0.) stop
00911 gamser=0.
00912 return
00913 endif
00914 ap=a
00915 sum=1./a
00916 del=sum
00917 do 11 n=1,itmax
00918 ap=ap+1.
00919 del=del*x/ap
00920 sum=sum+del
00921 if(abs(del).lt.abs(sum)*eps)go to 1
00922 11 continue
00923 write(0,*) 'aking: a too large, itmax too small'
00924 stop
00925 1 gamser=sum*exp(-x+a*log(x)-gln)
00926
00927 return
00928 end
00929
00930
00931 double precision function gammln(xx)
00932 implicit double precision (a-h,o-z)
00933 real*8 cof(6),stp,half,one,fpf,x,tmp,ser
00934 data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
00935 & -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
00936 data half,one,fpf/0.5d0,1.0d0,5.5d0/
00937
00938 x=xx-one
00939 tmp=x+fpf
00940 tmp=(x+half)*log(tmp)-tmp
00941 ser=one
00942 do 11 j=1,6
00943 x=x+one
00944 ser=ser+cof(j)/x
00945 11 continue
00946 gammln=tmp+log(stp*ser)
00947
00948 return
00949 end
00950
00951
00952 subroutine gcf(gammcf,a,x,gln)
00953 implicit double precision (a-h,o-z)
00954 parameter (itmax=100,eps=3.e-7)
00955
00956 gln=gammln(a)
00957 gold=0.
00958 a0=1.
00959 a1=x
00960 b0=0.
00961 b1=1.
00962 fac=1.
00963 do 11 n=1,itmax
00964 an=float(n)
00965 ana=an-a
00966 a0=(a1+a0*ana)*fac
00967 b0=(b1+b0*ana)*fac
00968 anf=an*fac
00969 a1=x*a0+anf*a1
00970 b1=x*b0+anf*b1
00971 if(a1.ne.0.)then
00972 fac=1./a1
00973 g=b1*fac
00974 if(abs((g-gold)/g).lt.eps)go to 1
00975 gold=g
00976 endif
00977 11 continue
00978 write(0,*) 'aking: a too large, itmax too small'
00979 stop
00980 1 gammcf=exp(-x+a*log(x)-gln)*g
00981
00982 return
00983 end
00984
00985
00986 double precision function notsogoodran2(idum)
00987
00988 c Notes from Steve, 7/16/98:
00989 c
00990 c This is the RAN2 random number generator from Numerical Recipes
00991 c (first FORTRAN edition, ~1986 or so). It has a number of limitations,
00992 c and seems to fail when too many calls are made to it..
00993 c
00994 c This generator will fail (with j out of range) if idum is too big.
00995 c For the parameters below, the limit is idum < 150999.
00996 c
00997 c (The mod function returns negative values for negative arguments
00998 c so, if ia*idum+ic <= 0 during the do-iteration below, j will go
00999 c out of range. This requires that the initial idum satisfy
01000 c
01001 c ia * (ic - idum) + ic > 0
01002 c or
01003 c idum < ic + ic/ia
01004 c = ic + 110
01005 c = 150999.)
01006
01007 implicit none
01008 integer m, ia, ic, iff, j, iy, idum
01009 integer idum_save, idum_prev, count
01010 double precision rm
01011 parameter (m=714025, ia=1366, ic=150889, rm=1.4005112e-6)
01012 integer ir(97)
01013
01014 c Note from Steve 7/15/98: the "save" here is crucially important!
01015
01016 save iff, ir, iy, idum_save, idum_prev, count
01017 data iff /0/ count /0/
01018
01019 idum_prev = idum_save
01020 idum_save = idum
01021 count = count + 1
01022
01023 if(idum.lt.0.or.iff.eq.0)then
01024 iff=1
01025 idum=mod(ic-idum, m)
01026 do j=1,97
01027 idum=mod(ia*idum+ic, m)
01028 ir(j)=idum
01029 enddo
01030 idum=mod(ia*idum+ic, m)
01031 iy=idum
01032 endif
01033 j=1+(97*iy)/m
01034 if(j.gt.97.or.j.lt.1) then
01035 write(0,*)'aking: ran2: j = ', j, ' count = ', count
01036 write(0,*)' idum = ', idum,
01037 + ' ia*idum+ic = ', ia*idum+ic
01038 write(0,*)' idum_prev = ', idum_prev,
01039 + ' ia*idum_prev+ic = ', ia*idum_prev+ic
01040 stop
01041 endif
01042 iy=ir(j)
01043
01044 notsogoodran2=iy*rm
01045
01046 idum=mod(ia*idum+ic, m)
01047 ir(j)=idum
01048
01049 c write(0,*)count, idum
01050
01051 return
01052 end
01053
01054
01055 double precision function betterran2(idum)
01056
01057 c Fortran-ized version of the ran2 function from "Numerical Recipes
01058 c in C," 2nd edition. Works much better than the older Fortran code,
01059 c and is also the basis of randinter() used in the rest of Starlab.
01060
01061
01062 c ***** An even better solution is to use a wrapper to allow *****
01063 c ***** Fortran to use the C++ function randinter! (Steve, 8/99) *****
01064
01065
01066 c Implemented by Steve, 8/99.
01067
01068 implicit none
01069
01070 integer idum
01071 integer IM1, IM2, IMM1, IA1, IA2, IQ1, IQ2, IR1, IR2, NTAB, NDIV
01072 double precision AM, EPS, RNMX
01073
01074 parameter (IM1 = 2147483563)
01075 parameter (IM2 = 2147483399)
01076 parameter (AM = (1.0/IM1))
01077 parameter (IMM1 = (IM1-1))
01078 parameter (IA1 = 40014)
01079 parameter (IA2 = 40692)
01080 parameter (IQ1 = 53668)
01081 parameter (IQ2 = 52774)
01082 parameter (IR1 = 12211)
01083 parameter (IR2 = 3791)
01084 parameter (NTAB = 32)
01085 parameter (NDIV = (1+IMM1/NTAB))
01086 parameter (EPS = 1.e-14)
01087 parameter (RNMX = (1.0-EPS))
01088
01089 integer init, j, k, idum2, iy, iv(NTAB)
01090
01091 save init, idum2, iy, iv
01092 data init /0/ idum2 /123456789/ iy /0/
01093
01094 c Initialize if idum < 0.
01095
01096 if (idum .le. 0 .or. init .eq. 0) then
01097
01098 c write(0,*)'Initializing....'
01099
01100 if (idum .lt. 0) idum = -idum
01101 if (idum .eq.0) idum = 1
01102 idum2 = idum
01103
01104 do j = NTAB+8, 1, -1
01105 k = (idum)/IQ1
01106 idum = IA1*(idum-k*IQ1) - k*IR1
01107 if (idum .lt. 0) idum = idum + IM1
01108 if (j .le. NTAB) iv(j) = idum
01109 end do
01110 iy = iv(1)
01111
01112 init = 1
01113 end if
01114
01115 c write(0,*)'idum = ', idum
01116
01117 k = (idum)/IQ1
01118 idum = IA1*(idum-k*IQ1) - k*IR1
01119 if (idum .lt. 0) idum = idum + IM1
01120
01121 c write(0,*)'idum = ', idum
01122
01123 k = idum2/IQ2
01124 idum2 = IA2*(idum2-k*IQ2)-k*IR2
01125 if (idum2 .lt. 0) idum2 = idum2 + IM2
01126
01127 c write(0,*)'idum2 = ', idum2
01128
01129 j = iy/NDIV + 1
01130 iy = iv(j) - idum2
01131 iv(j) = idum
01132 if (iy .lt. 1) iy = iy + IMM1
01133
01134 c write(0,*)'iy = ', iy
01135
01136 betterran2 = AM*iy
01137 if (betterran2 .gt. RNMX) betterran2 = RNMX
01138
01139 write(0,*) betterran2
01140
01141 return
01142 end
01143
01144
01145 subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,
01146 & derivs,rkqc)
01147 implicit double precision(a-h,o-z)
01148 parameter (maxstp=10000,nmax=10,two=2.0,zero=0.0,tiny=1.e-30)
01149 common /path/ kmax,kount,dxsav,xp(200),yp(10,200)
01150 dimension ystart(nvar),yscal(nmax),y(nmax),dydx(nmax)
01151 external derivs
01152
01153 x=x1
01154 h=sign(h1,x2-x1)
01155 nok=0
01156 nbad=0
01157 kount=0
01158 do 11 i=1,nvar
01159 y(i)=ystart(i)
01160 11 continue
01161 xsav=x-dxsav*two
01162 do 16 nstp=1,maxstp
01163 call derivs(x,y,dydx)
01164 do 12 i=1,nvar
01165 yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
01166 12 continue
01167 if(kmax.gt.0)then
01168 if(abs(x-xsav).gt.abs(dxsav)) then
01169 if(kount.lt.kmax-1)then
01170 kount=kount+1
01171 xp(kount)=x
01172 do 13 i=1,nvar
01173 yp(i,kount)=y(i)
01174 13 continue
01175 xsav=x
01176 endif
01177 endif
01178 endif
01179 if((x+h-x2)*(x+h-x1).gt.zero) h=x2-x
01180 call rkqc(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
01181 if(hdid.eq.h)then
01182 nok=nok+1
01183 else
01184 nbad=nbad+1
01185 endif
01186 if((x-x2)*(x2-x1).ge.zero)then
01187 do 14 i=1,nvar
01188 ystart(i)=y(i)
01189 14 continue
01190 if(kmax.ne.0)then
01191 kount=kount+1
01192 xp(kount)=x
01193 do 15 i=1,nvar
01194 yp(i,kount)=y(i)
01195 15 continue
01196 endif
01197 return
01198 endif
01199 if(abs(hnext).lt.hmin) then
01200 write(0,*) 'aking: stepsize smaller than minimum.'
01201 stop
01202 endif
01203 h=hnext
01204 16 continue
01205
01206 write(0,*) 'aking: too many steps.'
01207 stop
01208 end
01209
01210
01211 subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
01212 implicit double precision (a-h,o-z)
01213 parameter (nmax=10,fcor=.0666666667,
01214 & one=1.,safety=0.9,errcon=6.e-4)
01215 external derivs
01216 dimension y(n),dydx(n),yscal(n),ytemp(nmax),ysav(nmax),dysav(nmax)
01217
01218 pgrow=-0.20
01219 pshrnk=-0.25
01220 xsav=x
01221 do 11 i=1,n
01222 ysav(i)=y(i)
01223 dysav(i)=dydx(i)
01224 11 continue
01225 h=htry
01226 1 hh=0.5*h
01227 call rk4(ysav,dysav,n,xsav,hh,ytemp,derivs)
01228 x=xsav+hh
01229 call derivs(x,ytemp,dydx)
01230 call rk4(ytemp,dydx,n,x,hh,y,derivs)
01231 x=xsav+h
01232 if(x.eq.xsav) then
01233 write(0,*) 'aking: stepsize not significant in rkqc.'
01234 stop
01235 endif
01236 call rk4(ysav,dysav,n,xsav,h,ytemp,derivs)
01237 errmax=0.
01238 do 12 i=1,n
01239 ytemp(i)=y(i)-ytemp(i)
01240 errmax=max(errmax,abs(ytemp(i)/yscal(i)))
01241 12 continue
01242 errmax=errmax/eps
01243 if(errmax.gt.one) then
01244 h=safety*h*(errmax**pshrnk)
01245 goto 1
01246 else
01247 hdid=h
01248 if(errmax.gt.errcon)then
01249 hnext=safety*h*(errmax**pgrow)
01250 else
01251 hnext=4.*h
01252 endif
01253 endif
01254 do 13 i=1,n
01255 y(i)=y(i)+ytemp(i)*fcor
01256 13 continue
01257
01258 return
01259 end
01260
01261
01262 subroutine rk4(y,dydx,n,x,h,yout,derivs)
01263 implicit double precision (a-h,o-z)
01264 parameter (nmax=10)
01265 dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),dym(nmax)
01266
01267 hh=h*0.5
01268 h6=h/6.
01269 xh=x+hh
01270 do 11 i=1,n
01271 yt(i)=y(i)+hh*dydx(i)
01272 11 continue
01273 call derivs(xh,yt,dyt)
01274 do 12 i=1,n
01275 yt(i)=y(i)+hh*dyt(i)
01276 12 continue
01277 call derivs(xh,yt,dym)
01278 do 13 i=1,n
01279 yt(i)=y(i)+h*dym(i)
01280 dym(i)=dyt(i)+dym(i)
01281 13 continue
01282 call derivs(x+h,yt,dyt)
01283 do 14 i=1,n
01284 yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
01285 14 continue
01286
01287 return
01288 end
01289
01290
01291 subroutine derivs (x,y,dydx)
01292 implicit double precision (a-h,o-z)
01293 parameter (nvar=2,m=51)
01294 common /params/ w0,alpha1,alpha3,rho(m),rhodot(m),rlag,yking(2,m),
01295 & den,rking
01296 dimension y(nvar),dydx(nvar)
01297
01298 pi = 4.d0*datan(1.d0)
01299 if (x.ge.0.d0) then
01300 rhox =-sqrt(x)*(x+1.5d0)+0.75*sqrt(pi)*exp(x)
01301 & *dch_erf(sqrt(x))
01302 else
01303 rhox = 0.d0
01304 endif
01305 dydx(1) = y(2)
01306 dydx(2) = 0.25d0*y(2)**2*(6.d0+9.d0*y(2)*rhox/den)/y(1)
01307
01308 return
01309 end
01310
01311
01312 double precision function sigma(xstar,zstar,y)
01313 implicit double precision (a-h,o-z)
01314 parameter (ne=5,m=51)
01315 common /params/ w0,alpha1,alpha3,rho(m),rhodot(m),rlag,
01316 & yking(2,m),den,rking
01317 common x(m),h,mm,n,c2,anorm
01318 dimension y(ne,m)
01319
01320 densty(z) = -sqrt(z)*(z+1.5d0)+0.75*sqrt(pi)*exp(z)
01321 & *dch_erf(sqrt(z))
01322
01323 eps = -exp(y(5,m))
01324 s0 = 1.5d0*(1.d0+alpha3/alpha1)
01325 wc1 = -2.d0*yking(1,m)**1.5d0/yking(2,m)
01326 wc2 = -wc1/sqrt(yking(1,m))
01327 w0c1 = wc1*y(2,m)
01328 w0c2 = y(1,m) - w0c1/sqrt(yking(1,m))
01329 w2c1 = yking(1,m)**1.5d0*y(3,m)
01330 pi = 4.d0*atan(1.d0)
01331 sum = 0.d0
01332 ystar = 0.d0
01333 coeff = 0.5d0
01334 10 continue
01335 r2s2 = (alpha3/alpha1-0.5d0)*(zstar**2-0.5d0*(xstar**2+
01336 & ystar**2))+ 0.75d0*(xstar**2-ystar**2)
01337 r2 = xstar**2 + ystar**2 + zstar**2
01338 if (r2.gt.0.d0) then
01339 s2 = r2s2/r2
01340 else
01341 s2 = 1.d0
01342 endif
01343 r = sqrt(r2)
01344 if (r.lt.rking) then
01345 if (r.lt.sqrt(yking(1,1))) then
01346 w1 = x(1) + y(1,1)*s0 + y(3,1)*s2 - eps*r**2*
01347 & (2.d0/3.d0)*(s0/3.d0 + s2)
01348 w = w0 - (r2/yking(1,1))*(w0 - w1)
01349 else
01350 j = 1
01351 20 continue
01352 j = j + 1
01353 if (r.gt.sqrt(yking(1,j))) goto 20
01354 wj1 = x(j-1) + y(1,j-1)*s0 + y(3,j-1)*s2 - eps*r**2*
01355 & (2.d0/3.d0)*(s0/3.d0 + s2)
01356 wj = x(j) + y(1,j)*s0 + y(3,j)*s2 - eps*r**2*
01357 & (2.d0/3.d0)*(s0/3.d0 + s2)
01358 w = wj1 + (r2-yking(1,j-1))*(wj-wj1)/
01359 & (yking(1,j)-yking(1,j-1))
01360 endif
01361 else
01362 w = wc1/r + wc2
01363 & +s0*(w0c1/r + w0c2) + s2*w2c1/r**3 - eps*r**2*
01364 & (2.d0/3.d0)*(s0/3.d0 + s2)
01365 endif
01366 if (w.gt.0.d0) then
01367 rhost = densty(w)/den
01368 sum = sum + coeff*rhost
01369 if (coeff.lt.1.d0) coeff=1.d0
01370 ystar = ystar + 0.1
01371 goto 10
01372 else
01373 sigma = sum
01374 return
01375 endif
01376
01377 end
