program cjdet include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'network.dek' include 'cjdet.dek' c..evaluates a chapman-jouget detonation c..declare integer i,k,kkase double precision xh1,xhe4,xc12,xo16,xne22,qdum,mach_fac c..formats 01 format(1x,a) c..initialize the network call zet127 call init_torch c..get the upstream conditions 100 write(6,01) 'upstream temp, den, h1, he4, c12, o16 ne22 =>' read (5,*) temp_up,den_up,xh1,xhe4,xc12,xo16,xne22 c..see if we want a driven detonation write(6,*) ' ' write(6,*) 'to drive a detonation give the mach number' write(6,*) 'otherwise enter a number less than one =>' read(5,*) mach_fac if (mach_fac .lt. 1.0) mach_fac = 0.0d0 c..set the upsteam composition variables do i=1,ionmax xmass_up(i) = 1.0d-20 enddo if (ih1 .ne. 0) xmass_up(ih1) = xh1 if (ihe4 .ne. 0) xmass_up(ihe4) = xhe4 if (ic12 .ne. 0) xmass_up(ic12) = xc12 if (io16 .ne. 0) xmass_up(io16) = xo16 if (ine22 .ne. 0) xmass_up(ine22) = xne22 do i=1,ionmax xmass_up(i) = min(1.0d0,max(xmass_up(i),1.0d-20)) end do c..get the chapman-jouget detonation solution kkase = 1 mach = 0.0d0 call cjsolve(kkase,xmass_up,temp_up,den_up,mach, 1 qburn_cj,xmass_cj,ener_up,pres_up,cs_up, 2 vel_det,vel_cj,temp_cj,den_cj,ener_cj,pres_cj,cs_cj) c..get shock solution kkase = 4 mach_sh = vel_det/cs_up call cjsolve(kkase,xmass_up,temp_up,den_up,mach_sh, 1 qdum,xmass_up,ener_up,pres_up,cs_up, 2 vel_det,vel_sh,temp_sh,den_sh,ener_sh,pres_sh,cs_sh) mach = mach_fac * mach_sh if (mach .ne. 0.0) then c..get the shock values for this mach number kkase = 4 call cjsolve(kkase,xmass_up,temp_up,den_up,mach, 1 qdum,xmass_up,ener_up,pres_up,cs_up, 2 vel_det,vel_swsh,temp_swsh,den_swsh, 3 ener_swsh,pres_swsh,cs_swsh) c..get the strong point kkase = 2 call cjsolve(kkase,xmass_up,temp_up,den_up,mach, 1 qburn_s,xmass_s,ener_up,pres_up,cs_up, 2 vel_det,vel_s,temp_s,den_s,ener_s,pres_s,cs_s) c..get the weak point kkase = 3 call cjsolve(kkase,xmass_up,temp_up,den_up,mach, 1 qburn_w,xmass_w,ener_up,pres_up,cs_up, 2 vel_det,vel_w,temp_w,den_w,ener_w,pres_w,cs_w) end if c..write the out the results and go back for another call cjout(6) goto 100 end subroutine cjout(lu) include 'implno.dek' include 'vector_eos.dek' include 'network.dek' include 'cjdet.dek' c..writes the output to logical unit lu c..declare integer lu,i,indx(ionmax) double precision sum,dum,tau_qse,tau_nse,tau_burn, 1 abar,zbar c..formats 01 format(1x,a7,'=',1pe10.3,' ',a7,'=',1pe10.3,' ', 1 a7,'=',1pe10.3) 02 format(1x,a4,'=',1pe10.3,' ',a4,'=',1pe10.3,' ', 1 a4,'=',1pe10.3,' ',a4,'=',1pe10.3) 03 format(1x,a) write(lu,*) ' ' write(lu,03) 'upstream state:' write(lu,01) 'temp_up',temp_up,'den_up ',den_up, 1 'pres_up',pres_up write(lu,01) 'cs_up ',cs_up write(lu,*) ' ' write(lu,03) 'state behind shock:' write(lu,01) 'temp_sh',temp_sh,'den_sh ',den_sh, 1 'pres_sh',pres_sh write(lu,01) 'cs_sh ',cs_sh, 2 'vel_mat',vel_sh,'vel_det',vel_det write(lu,01) 'mach_sh',mach_sh write(lu,*) ' ' write(lu,03) 'cj state (should be sonic with vel_mat = cs_cj):' write(lu,01) 'temp_cj',temp_cj,'den_cj ',den_cj, 1 'pres_cj',pres_cj write(lu,01) 'cs_cj ',cs_cj, 2 'vel_mat',vel_cj,'vel_det',vel_det write(lu,01) 'mach_cj',vel_cj/cs_cj,'qburn_cj',qburn_cj c..nse composition at the cj point sum = 0.0d0 do i=1,ionmax sum = sum + xmass_cj(i) enddo sum = 1.0d0 - sum write(lu,*) ' ' write(lu,*) 'top 10 cj nse mass fractions:' call indexx(ionmax,xmass_cj,indx) write(lu,02) (ionam(indx(i)), 1 xmass_cj(indx(i)), i=ionmax,ionmax-9,-1), 2 '1-sum',sum c..qse and nse time and length scales tau_qse = exp(149.7d0/(temp_cj*1.0d-9) - 37.7d0) tau_nse = den_cj**(0.2d0)*exp(179.7d0/(temp_cj*1.0d-9) - 39.0d0) write(lu,*) ' ' write(lu,*) 'qse and nse time and length scales ' write(lu,01) 'tau_qse',tau_qse,'len_qse',tau_qse*vel_det write(lu,01) 'tau_nse',tau_nse,'len_nse',tau_nse*vel_det write(lu,*) ' ' c..the weak and strong detonations if (qburn_w .ne. 0.0) then write(lu,*) ' ' write(lu,*) 'strong and weak detonation shock state:' write(lu,01) 'temp_sw=',temp_swsh,' den_sw=',den_swsh, 1 ' pres_sw=',pres_swsh write(lu,01) 'cs_sw =',cs_swsh, 2 ' vel_sw=',vel_swsh,' vel_det=',vel_det write(lu,01) 'mach_sw=',mach c..strong detonations dum = 0.0d0 if (cs_s .ne. 0.0) dum = vel_s/cs_s write(lu,*) ' ' write(lu,*) 'strong detonation (subsonic with mach_s < 1):' write(lu,01) 'temp_s =',temp_s,' den_s =',den_s, 1 ' pres_s =',pres_s write(lu,01) 'cs_s =',cs_s, 2 ' vel_s =',vel_s,' vel_det=',vel_det write(lu,01) 'mach_s =',dum,' qburn_s =',qburn_s sum = 0.0d0 do i=1,ionmax sum = sum + xmass_s(i) enddo sum = 1.0d0 - sum write(lu,*) ' ' write(lu,*) 'top 10 strong detonation nse mass fractions:' call indexx(ionmax,xmass_s,indx) write(lu,02) (ionam(indx(i)), 1 xmass_s(indx(i)), i=ionmax,ionmax-9,-1), 2 '1-sum',sum write(lu,*) ' ' c..weak detonations dum = 0.0d0 if (cs_w .ne. 0.0) dum = vel_w/cs_w write(lu,*) ' ' write(lu,*) 'weak detonation (supersonic mach_w > 1):' write(lu,01) 'temp_w =',temp_w,' den_w =',den_w, 1 ' pres_w =',pres_w write(lu,01) 'cs_w =',cs_w, 2 ' vel_w =',vel_w,' vel_det=',vel_det write(lu,01) 'mach_w =',dum,' qburn_w =',qburn_w sum = 0.0d0 do i=1,ionmax sum = sum + xmass_w(i) enddo sum = 1.0d0 - sum write(lu,*) ' ' write(lu,*) 'top 10 weak detonation nse mass fractions:' call indexx(ionmax,xmass_w,indx) write(lu,02) (ionam(indx(i)), 1 xmass_w(indx(i)), i=ionmax,ionmax-9,-1), 2 '1-sum',sum write(lu,*) ' ' end if return end subroutine cjsolve(kkase,xmass_up,temp_up,den_up,mach, 1 qx,xmass_det,ener_up,pres_up,cs_up, 2 vel_det,vel_x,temp_x,den_x,ener_x,pres_x,cs_x) include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'network.dek' c..solves the hugoniot and rayleigh relations for a detonation or a shock. c..an nse distribution is assumed for the detonated material. c..input: c..kkase = 1 for a chapman-jouget detonation c.. = 2 for a strong point driven detonation c.. = 3 for a weak point driven detonation c.. = 4 for a shock wave c..xmass_up = upstream composition vector c..temp_up = temperature of upstream material c..den_up = density of upstream material c..mach = mach number of shock or detonation c..output: c..qx = energy/gram from burning c..xmass_det = burned composition c..ener_up = energy of upstream material c..pres_up = pressure of upstream material c..cs_up = sound speed of upstream material c..vel_det = speed of detonation or shock front = upstream speed c..vel_x = speed behind detonation or shock = downstream speed c..temp_x = temperature of downstream material c..den_x = density of downstream material c..ener_x = internal energy of downstream material c..pres_x = pressure of downstream material c..cs_x = sound speed of downstream material c..declare the pass integer kkase double precision xmass_up(1),temp_up,den_up,mach,qx,xmass_det(1), 1 ener_up,pres_up,cs_up,vel_det,vel_x,temp_x,den_x, 2 ener_x,pres_x,cs_x c..common block communication with the routine cjfunc integer kase double precision xmup(abignet),den1,temp1,p1,u1,v1,ye1,vel1,cs1, 1 mach1,qburn,den2,temp2,p2,u2,v2,vel2,cs2 common /cjstate/ xmup,den1,temp1,p1,u1,v1,ye1,vel1,cs1, 1 mach1,qburn,den2,temp2,p2,u2,v2,vel2,cs2, 2 kase c..locals external cjfunc character*8 type logical check integer i,n,ntrial,ntaken,nfev,nstrong,nsmax,nweak,nwmax parameter (ntrial = 60, n=2, nsmax = 10, nwmax=10) double precision cv1,g1,g1p1,g1m1,msq,v1mv2,x(n),f(n),dum, 1 xl,xx,gsq,den2_sav,abar,zbar,wbar,ye,xcess, 2 tolx,tolf,conv parameter (tolf = 1.0d-6, tolx = 1.0d-6) parameter (conv = ev2erg * 1.0d6 * avo) c..check the input if (kkase .lt. 1 .or. kkase .gt. 4) then write(6,*) 'kkase =',kkase stop 'kkase in cjsolve is invalid' end if c..transfer the input to common kase = kkase temp1 = temp_up den1 = den_up mach1 = mach do i=1,ionmax xmup(i) = xmass_up(i) enddo c..load the eos conditions call azbar(xmass_up,aion,zion,wion,ionmax, 1 ymass,abar,zbar,wbar,ye,xcess) temp_row(1) = temp_up den_row(1) = den_up abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 c..call an eos call helmeos c..set upstream thermodynamic conditions p1 = ptot_row(1) u1 = etot_row(1) cs1 = cs_row(1) cv1 = cv_row(1) g1 = gam1_row(1) v1 = 1.0d0/den1 ye1 = zbar/abar pres_up = p1 ener_up = u1 cs_up = cs1 c..now we start making initial guesses for downstream temperature and density c..for the detonation cases here is a guess for the energy generated c..a pure si28 burned composition seems robust if (kase .eq. 1 .or. kase .eq. 2 .or. kase .eq. 3) then qburn = (1.0d0 - xmass_up(isi28))/aion(isi28)*bion(isi28)*conv end if c..an initial guess for the density and temperature of the burned material c..from landau & lifshitz fluid dynamics, 129.15 c..with gamma2 about gamma1 and cv2 about 1/8 cv1 c..this is nearly exact for the ions if (kase .eq. 1 .or. kase .eq. 2 .or. kase .eq. 3) then g1p1 = g1 + 1.0d0 den2 = den1 * g1p1/g1 temp2 = 2.0d0 * g1 * qburn / (8.0d0 * cv1 *g1p1) c..limit temp2 so nse material is not photodisintegrated back to helium temp2 = min(5.0d9, max(temp1,temp2)) c..modify density guess for strong or weak detonations c..the weak detonation appears to be a much stronger attractor if (kase .eq. 2) then den2 = 1.2d0 * den2 else if (kase .eq. 3) then den2 = 0.9 * den2 end if c..an initial guess for the density and temperature behind a shock wave c..from landau & lifshitz fluid dynamics, eq 85.7 to 85.10 c..its exact for the ions else g1p1 = g1 + 1.0d0 g1m1 = g1 - 1.0d0 msq = mach1*mach1 den2 = den1 * g1p1*msq/(g1m1*msq + 2.0d0) temp2 = temp1 * (2.0d0*g1*msq - g1m1) * 1 (g1m1*msq + 2.0d0) / (g1p1*g1p1*msq) end if c..loop to here with a new den2 if the strong or weak point is not proper nstrong = 0 nweak = 0 den2_sav = den2 100 continue c..done with the initial guesses section c..root find of the rayleigh line and the hugoniot x(1) = den2 x(2) = temp2 call xnewt_cj(ntrial,x,n,tolx,tolf,ntaken,check,nfev,cjfunc) if (check .or. ntaken .eq. ntrial) then write(6,*) write(6,*) 'check convergence of xnewt_cj root find' write(6,*) write(6,*) 'iterations taken =',ntaken write(6,*) 'function evals =',nfev write(6,111) 'roots =',x(1),x(2) 111 format(1x,a,1p2e14.6) end if c..with the converged values, get the return arguments call cjfunc(dum,x,f) c..set the return arguments do i=1,ionmax xmass_det(i) = xmass(i) enddo vel1 = den2/den1 * vel2 vel_det = vel1 vel_x = vel2 temp_x = temp2 den_x = den2 pres_x = p2 ener_x = u2 cs_x = cs2 qx = qburn c..the strong point solution needs to be checked before returning if (kase .eq. 2 .and. nstrong .gt. nsmax) then write(6,*) ' ' write(6,*) 'warning: did not find strong point solution' write(6,*) ' ' return end if if (kase .eq. 2) then if (vel_x .ge. cs_x) then nstrong = nstrong + 1 den2 = 1.1d0 * den2_sav den2_sav = den2 goto 100 end if end if c..the weak point solution needs to be checked before returning if (kase .eq. 3 .and. nweak .gt. nwmax) then write(6,*) ' ' write(6,*) 'warning: did not find weak point solution' write(6,*) ' ' return end if if (kase .eq. 3) then if (vel_x .le. cs_x) then nweak = nweak + 1 den2 = 0.9d0 * den2_sav den2_sav = den2 goto 100 end if end if c..normal bail point return end subroutine cjfunc(x,y,f) include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'network.dek' c..this routine returns the functions to do the root find on c..input is x (not relevant here) and y, a vector of the unknowns, c..y(1) is the density and y(2) is the temperature. c..output is f, a vector of residuals to be minimized. f(1) is the c..rayleigh line and f(2) is the hugoniot. c..declare the pass double precision x,y(*),f(*) c..common block communication with the routine cjfunc integer kase double precision xmup(abignet),den1,temp1,p1,u1,v1,ye1,vel1,cs1, 1 mach1,qburn,den2,temp2,p2,u2,v2,vel2,cs2 common /cjstate/ xmup,den1,temp1,p1,u1,v1,ye1,vel1,cs1, 1 mach1,qburn,den2,temp2,p2,u2,v2,vel2,cs2, 2 kase c..locals integer i,newguess,iprint double precision xmnse(abignet),xmunn,xmupp,abar,zbar,wbar, 1 ye,xcess,conv parameter (conv = ev2erg * 1.0d6 * avo) c..map the input vector, bail if its hosed den2 = y(1) temp2 = y(2) if (den2 .lt. 0.0 .or. temp2 .le. 0.0) then f(1) = 1.0d30 f(2) = 1.0d30 return end if c..load the nse composition if (kase .eq. 1 .or. kase .eq. 2 .or. kase .eq. 3) then iprint = 0 newguess = 1 call nse(temp2,den2,ye1,newguess,1,1,xmnse,xmunn,xmupp,iprint) end if c..get the energy generated qburn = 0.0d0 if (kase .eq. 1 .or. kase .eq. 2 .or. kase .eq. 3) then do i=1,ionmax qburn = qburn + (xmnse(i) - xmup(i))/aion(i) *bion(i) end do qburn = qburn * conv end if c..get the eos c..for the detonation cases use the nse composition c..for the shock case, the shocked composition is the upstream composition if (kase .eq. 1 .or. kase .eq. 2 .or. kase .eq. 3) then do i=1,ionmax xmass(i) = xmnse(i) enddo else do i=1,ionmax xmass(i) = xmup(i) enddo end if call azbar(xmass,aion,zion,wion,ionmax, 1 ymass,abar,zbar,wbar,ye,xcess) temp_row(1) = temp2 den_row(1) = den2 abar_row(1) = abar zbar_row(1) = zbar jlo_eos = 1 jhi_eos = 1 call helmeos c..set the downstream thermodynamic variables v2 = 1.0d0/den_row(1) p2 = ptot_row(1) u2 = etot_row(1) cs2 = cs_row(1) c..for a chapman-jouget detonation vel2 is the burned sound speed c..for a strong or weak detonation or a shock wave the upstream c..mach number mach1 is specified if (kase .eq. 1) then vel2 = cs2 else vel2 = mach1 * cs1 * den1/den2 endif c..the rayleigh line, glassman, page 227, eq. 6, fickett & davis page 17 f(1) = (p2 - p1) - (vel2*vel2*den2*den2) * (v1 - v2) c..the specific internal energy hugoniot, glassman, page 229, eq. 11 f(2) = u1 + qburn - u2 + 0.5d0 *( (p1+p2) * (v1-v2) ) c..scale the functions for better behavior of the root finder f(1) = f(1) * v1/p1 f(2) = f(2) / (p1 * v1) return end subroutine xnewt_cj(ntrial,x,n,tolx,tolf,ntaken,check,nfev,func) include 'implno.dek' c..given an initial guess x(1:n) for the root of n equations, this routine c..finds the root by a globally convergent newtons method. the vector of c..functions to be zeroed, called fvec(1:n) in the routine below, is c..returned by the user supplied routine func. the output quantity check c..is false on nomal return and true if the routine has converged to a c..local minimum of the function xfminx_cj. if so, try restarting from a c..different initial guess. c.. c..np is the maximum number of equations n c..ntrial is the maximum number of iterations to try c..ntaken is the number of iterations done c..tolf sets the convergence on function values c..tolmin sets the criterion for deciding wheather spurious convergence to c.. a false minimum of xfminx_cj has occured c..tolx is the convergence criteria on deltax c..stpmx is the scaled maximum step length allowed in the line searches c..nfev is the number of function evaluations c..declare the pass external func logical check integer ntrial,n,ntaken,nfev double precision x(n),tolf,tolx c..locals integer np parameter (np=4) integer nn,i,its,j,indx(np) double precision fvec(np),tolmin,stpmx,d,den,f, 1 fold,stpmax,sum,temp,test,fjac(np,np),g(np), 2 p(np),xold(np),xfminx_cj,dum parameter (tolmin = 1.0d-12, 1 stpmx = 2.0d0) c..common block communicates values from routine xfminx_cj common /newtcj/ fvec,nn c..initialize if (n .gt. np) stop 'n > np in routine xnewt' nn = n f = xfminx_cj(x,func) nfev = 1 c.. test for the initial guess being a root, using a more stringent tolf test = 0.0d0 do i=1,n if (abs(fvec(i)) .gt. test) test = abs(fvec(i)) enddo if (test .lt. 0.01*tolf) then check = .false. return end if c..get stpmax for the line search sum = 0.0d0 do i=1,n sum = sum + x(i)*x(i) enddo stpmax = stpmx * max(sqrt(sum),dfloat(n)) c..start of iteration loop; get the jacobian do its = 1, ntrial ntaken = its c..second order accurate jacobian call jac_cj(dum,x,fjac,n,n,np,np,func) nfev = nfev + 2*n + 1 c..compute grad f for the line searches do i=1,n sum = 0.0d0 do j=1,n sum = sum + fjac(j,i)*fvec(j) enddo g(i) = sum enddo c..store x, and f and form right hand sides do i=1,n xold(i) = x(i) enddo fold = f do i=1,n p(i) = -fvec(i) enddo c..solve the linear systems call ludcmp_cj(fjac,n,np,indx,d) call lubksb_cj(fjac,n,np,indx,p) c..line search returns new x and f c..it also gets fvec at the new x when it calls xfminx_cj call lnsrch_cj(n,xold,fold,g,p,x,f,stpmax,check,nfev,func) c..test for convergence on function value test = 0.0d0 do i=1,n if (abs(fvec(i)) .gt. test) test = abs(fvec(i)) enddo if (test .lt. tolf) then check = .false. return end if c..check for zero gradiant, i.e. spurious convergence if (check) then test = 0.0d0 den = max(f, 0.5d0 * n) do i=1,n temp = abs(g(i)) * max(abs(x(i)),1.0d0)/den if (temp .gt. test) test = temp enddo if (test .lt. tolmin) then check = .true. else check = .false. end if return end if c..test for convergence on deltax test = 0.0d0 do i=1,n temp = (abs(x(i)-xold(i)))/max(abs(x(i)),1.0d0) if (temp .gt. test) test = temp enddo if (test .lt. tolx) return c..back for another iteration enddo check = .true. return end subroutine lnsrch_cj(n,xold,fold,g,p,x,f,stpmax,check,nfev,func) include 'implno.dek' c..given an n dimensional point xold(1:n), the value of the function fold c..and the gradient g(1:n) at the point, and a direction p(1:n), this routine c..finds a new point x(1:n) along the direction of p from xold where the c..function xfminx_cj has decreased "sufficiently". the new function value is c..returned in f. stpmax is an input quanity that limits the length of the c..steps so that the function is not evaluated in regions where it is c..undefined or subject to overflow. p is usually the newton direction. the c..output quantity check is false on normal exit, and true when x is too c..close to xold. in a minimization routine, this usually signals c..convergence and can be ignored. however, in a root finding routine, the c..calling routine should check wheather the convergence is spurious. c..declare the pass external func logical check integer n,nfev double precision f,fold,stpmax,g(n),p(n),x(n),xold(n) c..locals integer i double precision xfminx_cj,alf,tolx,a,alam,alam2,alamin,b, 1 disc,f2,rhs1,rhs2,slope,sum,temp,test,tmplam parameter (alf=1.0d-4, tolx=3.0d-13) c..alf ensures sufficient decrease in the function value c..tolx is the convergence criterion on deltax c..initialize and scale if the attempted step is too big check = .false. sum = 0.0d0 do i=1,n sum = sum + p(i)*p(i) enddo sum = sqrt(sum) if (sum .gt. stpmax) then do i=1,n p(i) = p(i) * stpmax/sum enddo end if slope = 0.0d0 do i=1,n slope = slope + g(i)*p(i) enddo if (slope .ge. 0.0) stop 'roundoff problem in lnsrch_cj' c..compute lambda_min test = 0.0d0 do i=1,n temp = abs(p(i))/max(abs(xold(i)),1.0d0) if (temp .gt. test) test = temp enddo alamin = tolx/test c..always try a full newton step, start of iteration loop alam = 1.0d0 1 continue do i=1,n x(i) = xold(i) + alam*p(i) enddo f = xfminx_cj(x,func) nfev = nfev + 1 c..convergence on deltax, for root finding, the calling routine c..should verify the convergence if (alam .lt. alamin) then do i=1,n x(i) = xold(i) enddo check = .true. return c..sufficient function decrease else if (f .le. fold + alf*alam*slope) then return c..backtrack else if (alam .eq. 1.0) then tmplam = -slope / (2.0d0 * (f-fold-slope)) else rhs1 = f - fold - alam*slope rhs2 = f2 - fold - alam2*slope a = (rhs1/alam**2 - rhs2/alam2**2)/(alam-alam2) b = (-alam2*rhs1/alam**2 + alam*rhs2/alam2**2) / (alam-alam2) if (a .eq. 0.0) then tmplam = -slope/(2.0d0 * b) else disc = b*b - 3.0d0 * a * slope if (disc .lt. 0.0) then tmplam = 0.5d0 * alam else if (b .le. 0.0) then tmplam = (-b + sqrt(disc)) / (3.0d0 * a) else tmplam = -slope/(b + sqrt(disc)) end if end if if (tmplam .gt. 0.5d0*alam) tmplam = 0.5d0*alam end if end if c..store for the next trip through alam2 = alam f2 = f alam = max(tmplam, 0.1d0*alam) go to 1 end double precision function xfminx_cj(x,func) include 'implno.dek' c..returns f = 0.5 f dot f at x. func is a user supplied routine of the c..functions to be root found. c.. c..common block communicates values back to routine xnewt c..declare external func integer nn,np,i parameter (np = 4) double precision x(1),fvec(np),sum,dum common /newtcj/ fvec,nn call func(dum,x,fvec) sum = 0.0d0 do i=1,nn sum = sum + fvec(i)*fvec(i) enddo xfminx_cj = 0.5d0 * sum return end subroutine jac_cj(x,y,dfdy,mcol,nrow,mmax,nmax,derivs) include 'implno.dek' c..this routine computes a second order accurate jacobian matrix c..of the function contained in the routine derivs. c.. c..input is the point x and the the vector y at which to compute the c..jacobian dfdy. c.. c..y has logical dimension nrow and physical dimension nmax, c..dfdy has logical dimension (mcol,nrow) and physical dimension (mmax,nmax) c.. c..uses 2*nrow + 1 function evaluations c..declare the pass external derivs integer mcol,nrow,mmax,nmax double precision x,y(nmax),dfdy(mmax,nmax) c..locals integer i,j,imax parameter (imax = 4) double precision fminus(imax),fplus(imax),rel,ax,temp,h,hinv parameter (rel = 3.162278d-8, ax = 1.0e-16) c..check if (nrow .gt. imax) stop 'nrow > imax in jac_cj' c..for each row, get the right stepsize do j=1,nrow temp = y(j) h = rel * max(abs(y(j)),ax) y(j) = temp + h h = y(j) - temp call derivs(x,y,fplus) y(j) = temp temp = y(j) y(j) = temp - h h = temp - y(j) call derivs(x,y,fminus) y(j) = temp c..compute the jth row of the jacobian hinv = 1.0d0/(2.0d0 * h) do i=1,mcol dfdy(i,j) = (fplus(i) - fminus(i)) * hinv enddo enddo c..restore the original state call derivs(x,y,fplus) return end subroutine zet127 include 'implno.dek' include 'network.dek' c.. c..this routine sets up a 127 isotope torch network c.. c..set the network id and name idnet = idtorch127 netname = 'torch127' inetin = 0 c..deuterium and tritium inetin = inetin + 1 izzz(inetin) = 1 inmin(inetin) = 1 inmax(inetin) = 2 c..helium 3 inetin = inetin + 1 izzz(inetin) = 2 inmin(inetin) = 1 inmax(inetin) = 1 c..lithium 7 inetin = inetin + 1 izzz(inetin) = 3 inmin(inetin) = 4 inmax(inetin) = 4 c..berylium 7-8 inetin = inetin + 1 izzz(inetin) = 4 inmin(inetin) = 3 inmax(inetin) = 3 c..boron 8 inetin = inetin + 1 izzz(inetin) = 5 inmin(inetin) = 3 inmax(inetin) = 3 c..carbon 12-14 inetin = inetin + 1 izzz(inetin) = 6 inmin(inetin) = 6 inmax(inetin) = 8 c..nitrogen 13-15 inetin = inetin + 1 izzz(inetin) = 7 inmin(inetin) = 6 inmax(inetin) = 8 c..oxygen 14-19 inetin = inetin + 1 izzz(inetin) = 8 inmin(inetin) = 6 inmax(inetin) = 11 c..flourine 17-21 inetin = inetin + 1 izzz(inetin) = 9 inmin(inetin) = 8 inmax(inetin) = 12 c..neon 18-24 inetin = inetin + 1 izzz(inetin) = 10 inmin(inetin) = 8 inmax(inetin) = 14 c..sodium 19-24 inetin = inetin + 1 izzz(inetin) = 11 inmin(inetin) = 8 inmax(inetin) = 13 c..magnesium 22-27 inetin = inetin + 1 izzz(inetin) = 12 inmin(inetin) = 10 inmax(inetin) = 15 c..aluminum 25-29 inetin = inetin + 1 izzz(inetin) = 13 inmin(inetin) = 12 inmax(inetin) = 16 c..silicon 27-30 inetin = inetin + 1 izzz(inetin) = 14 inmin(inetin) = 13 inmax(inetin) = 16 c..phosphorus 29-32 inetin = inetin + 1 izzz(inetin) = 15 inmin(inetin) = 14 inmax(inetin) = 18 c..sulfer 30-34 inetin = inetin + 1 izzz(inetin) = 16 inmin(inetin) = 14 inmax(inetin) = 18 c..clorine 33-37 inetin = inetin + 1 izzz(inetin) = 17 inmin(inetin) = 16 inmax(inetin) = 20 c..argon 35-39 inetin = inetin + 1 izzz(inetin) = 18 inmin(inetin) = 17 inmax(inetin) = 21 c..potasium 37-42 inetin = inetin + 1 izzz(inetin) = 19 inmin(inetin) = 18 inmax(inetin) = 23 c..calcium 39-44 inetin = inetin + 1 izzz(inetin) = 20 inmin(inetin) = 19 inmax(inetin) = 24 c..scandium 42-46 inetin = inetin + 1 izzz(inetin) = 21 inmin(inetin) = 21 inmax(inetin) = 25 c..titanium 43-49 inetin = inetin + 1 izzz(inetin) = 22 inmin(inetin) = 21 inmax(inetin) = 27 c..vandium 46-51 inetin = inetin + 1 izzz(inetin) = 23 inmin(inetin) = 23 inmax(inetin) = 28 c..chromium 47-52 inetin = inetin + 1 izzz(inetin) = 24 inmin(inetin) = 23 inmax(inetin) = 28 c..manganese 49-55 inetin = inetin + 1 izzz(inetin) = 25 inmin(inetin) = 24 inmax(inetin) = 30 c..iron 51-56 inetin = inetin + 1 izzz(inetin) = 26 inmin(inetin) = 25 inmax(inetin) = 30 c..cobolt 53-58 inetin = inetin + 1 izzz(inetin) = 27 inmin(inetin) = 26 inmax(inetin) = 31 c..nickel 54-60 inetin = inetin + 1 izzz(inetin) = 28 inmin(inetin) = 26 inmax(inetin) = 32 return end subroutine init_torch include 'implno.dek' include 'const.dek' include 'network.dek' c.. c..this routine initializes stuff for the torch network c.. c..declare logical ibhere integer zmax,fil14 parameter (zmax=85) parameter (fil14=14) character*132 string,word character*2 zsymb(zmax) character*10 aname integer i,j,k,l,ii,jj,ll,llmin,llmax,mm,j1,nz,na,nn,kk, 1 aidmin(zmax),aidmax(zmax),inta,intz double precision mev2erg,mev2gr parameter (mev2erg = ev2erg*1.0d6, 1 mev2gr = mev2erg/clight**2) c..for easy zeroing of the isotope pointers integer isotp(nisotp) equivalence (isotp(1),ih1) c..for easy zeroing of the rate pointers integer rts(numrates) equivalence (rts(1),ir3a) double precision dum,qful,xx c..here are the root isotope names data zsymb/'h ','he','li','be','b ','c ','n ','o ','f ','ne', 1 'na','mg','al','si','p ','s ','cl','ar','k ','ca', 2 'sc','ti','v ','cr','mn','fe','co','ni','cu','zn', 3 'ga','ge','as','se','br','kr','rb','sr','y ','zr', 4 'nb','mo','tc','ru','rh','pd','ag','cd','in','sn', 5 'sb','te','i' ,'xe','cs','ba','la','ce','pr','nd', 6 'pm','sm','eu','gd','tb','dy','ho','er','tm','yb', 7 'lu','hf','ta','w' ,'re','os','ir','pt','au','hg', 8 'tl','pb','bi','po','at'/ c..here are the min and max a's for each z data aidmin/ 2, 3, 6, 7, 8, 9, 11, 13, 14, 16, 1 17, 18, 20, 22, 23, 24, 25, 27, 30, 30, 2 34, 34, 38, 38, 42, 42, 46, 46, 50, 51, 3 55, 55, 59, 59, 63, 63, 68, 68, 72, 72, 4 76, 77, 81, 81, 85, 86, 88, 90, 92, 94, 5 97, 99, 101, 103, 106, 108, 110, 113, 115, 118, 6 120, 123, 125, 128, 130, 133, 136, 138, 141, 143, 7 146, 150, 153, 154, 160, 160, 164, 165, 168, 170, 8 172, 174, 176, 182, 188/ data aidmax/ 3, 6, 9, 12, 14, 18, 21, 22, 26, 31, 1 44, 47, 51, 54, 57, 60, 63, 67, 70, 73, 2 76, 80, 83, 86, 89, 92, 96, 99, 102, 105, 3 108, 112, 115, 118, 121, 124, 128, 131, 134, 137, 4 140, 144, 147, 150, 153, 156, 160, 163, 166, 169, 5 171, 173, 175, 177, 179, 181, 183, 185, 187, 189, 6 191, 193, 195, 197, 199, 201, 203, 205, 207, 209, 7 213, 214, 219, 220, 225, 226, 231, 234, 237, 240, 8 245, 246, 251, 237, 239/ c..popular format statements 01 format(a,i4) 06 format(2i5,f10.4) 07 format(f10.3) 08 format(6f10.3) 09 format(1x,i4,i4,i4,' ',a5) c..zero all the isotope pointers do i=1,nisotp isotp(i) = 0 enddo c..zero all the rate pointers do i=1,numrates rts(i) = 0 enddo c..decide on the arrow orientation c..downarrow true puts neutron, protons, alfa to the end c..downarrow false (i.e uparrow) puts neutron, protons, alfa at the beginning c..in general downarrowtrue is faster for dense linear algebra, and either c..orientation for sparse linear algebra. gift routines, however, do c..much better with uparrow (downarrow false). downarrow = .true. c downarrow = .false. c..set the beginning isotope index if (downarrow) then ionbeg = 1 else ionbeg = 4 end if c..open the nuclear recation rate data file c..use a soft link to connect bdat to the desired burn data file open(unit=fil14,file='BDAT',status='old') c..now start reading the nuclear reaction rate data file c..i is the code number of element z(i),n(i). c..j = 1 = ng j = 6 = an c..j = 2 = pn j = 7 = ag c..j = 3 = ground state b- j = 8 for semi-empirical electron captur c..j = 4 = pg j = 9 for semi-empirical positron decay c..j = 5 = ap j = 10 for semi-empirical beta decay c..ic1(j,i) = type formula to be used to calculate rate c..ic2(j,i) = number of constants in fitting reaction j on species i c..ic3(j,i) = where to start counting ic2 from c..initialize counters nful = 0 nfulnot = 0 k = 1 i = ionbeg - 1 do nn=1,ionmax icode2(nn) = 0 enddo c..put neutrons, protons and alfa first for up-arrow orientation if (.not.downarrow) then ineut = 1 aion(ineut) = 1.0d0 nion(ineut) = 1.0d0 zion(ineut) = 0.0d0 bion(ineut) = 0.0d0 c mion(ineut) = nion(ineut)*mn + zion(ineut)*mp + zion(ineut)*me c 1 -bion(ineut)*mev2gr mion(ineut) = nion(ineut)*mn +zion(ineut)*mp -bion(ineut)*mev2gr wion(ineut) = avo * mion(ineut) wpart(ineut) = 2.0d0 ionam(ineut) = 'neut' iprot = 2 aion(iprot) = 1.0d0 nion(iprot) = 0.0d0 zion(iprot) = 1.0d0 bion(iprot) = 0.0d0 c mion(iprot) = nion(iprot)*mn + zion(iprot)*mp + zion(iprot)*me c 1 - bion(iprot)*mev2gr mion(iprot) = nion(iprot)*mn +zion(iprot)*mp -bion(iprot)*mev2gr wion(iprot) = avo * mion(iprot) wpart(iprot) = 2.0d0 ionam(iprot) = 'prot' ihe4 = 3 aion(ihe4) = 4.0d0 nion(ihe4) = 2.0d0 zion(ihe4) = 2.0d0 bion(ihe4) = 28.29603d0 c mion(ihe4) = nion(ihe4)*mn + zion(ihe4)*mp + zion(ihe4)*me c 1 - bion(ihe4)*mev2gr mion(ihe4) = nion(ihe4)*mn + zion(ihe4)*mp - bion(ihe4)*mev2gr wion(ihe4) = avo * mion(ihe4) wpart(ihe4) = 1.0d0 ionam(ihe4) = 'he4' end if c..we keep returning here from various goto and loop constructions 60 i = i+1 c..read in the z and a and any fitting constants c write(6,*) 'reading', i read(fil14,02) nz,na,(ic1(j,i),ic2(j,i), j=1,10) 02 format(2i6,20i3) c write(6,*) 'read', nz,na zion(i) = nz aion(i) = na if (nz .eq. 99) go to 120 c..temperature dependent partition function information llmin = 5*(i-1)+1 llmax = llmin + 4 if (llmax .gt. 6*abignet) stop 'past as bounds' read(fil14,03) nz,nn,bion(i),(as(ll),ll=llmin,llmax),ist(i),aname 03 format(2i3,f11.4,f5.1,4e12.3,i2,a10) nion(i) = nn if (ist(i).ne.0) then if (6*i-6+2*ist(i) .gt. 6*abignet) stop 'past gs bounds' read(fil14,04) (gs(ll),ll=6*i-5,6*i-6+2*ist(i)) end if 04 format(f10.4,f10.3,f10.4,f10.3,f10.4,f10.3,f10.4,f10.3) c..decide if this isotope is in the network and branch accordingly do jj=1,inetin if (int(zion(i)) .eq. izzz(jj) .and. 1 int(nion(i)) .ge. inmin(jj) .and. 2 int(nion(i)) .le. inmax(jj)) goto 90 enddo c..not using this isotope, but c..do the read, backup i by one, and go back to 60 for another isotope do jj=1,10 if (ic1(jj,i) .gt. 0) then c..bdat921 takes format 05, rath_005.bdat takes format 240 c read(fil14,05) (dum,j1=1,ic2(jj,i)) read(fil14,240) (dum,j1=1,ic2(jj,i)) 05 format(7e10.3) 240 format (7e13.6) end if enddo i = i - 1 go to 60 c..using this isotope, read parameters for reaction j on species i 90 continue if (i .gt. abignet) stop 'abignet too small in init_torch' c..here are the isotopes we are using c call sqeeze(aname) c write(6,117) i,int(zion(i)),int(aion(i)),aname c 117 format(1x,3i4,' ',a) do j=1,8 ic3(j,i) = k if (ic1(j,i) .gt. 0) then kmax = k + ic2(j,i)-1 if (kmax .gt. cxdim) then write(6,*) 'kmax =', kmax,' cxdim =',cxdim stop 'kmax > cxdim in routine init_torch' end if c..bdat921 takes format 05, rath_005.bdat takes format 240 c read(fil14,05) (cx(j1), j1=k,kmax) read(fil14,240) (cx(j1), j1=k,kmax) k = kmax + 1 end if enddo c..and go back for another isotope go to 60 c..all done with the this part of the loading 120 continue c..append neutrons, protons and alfa if down arrow c..set the ending isotope index ionend if (downarrow) then ionmax = i ineut = ionmax aion(ineut) = 1.0d0 nion(ineut) = 1.0d0 zion(ineut) = 0.0d0 bion(ineut) = 0.0d0 c mion(ineut) = nion(ineut)*mn + zion(ineut)*mp + zion(ineut)*me c 1 - bion(ineut)*mev2gr mion(ineut) = nion(ineut)*mn +zion(ineut)*mp -bion(ineut)*mev2gr wion(ineut) = avo * mion(ineut) wpart(ineut) = 2.0d0 ionam(ineut) = 'neut' c.. ionmax = ionmax + 1 iprot = ionmax aion(iprot) = 1.0d0 nion(iprot) = 0.0d0 zion(iprot) = 1.0d0 bion(iprot) = 0.0d0 c mion(iprot) = nion(iprot)*mn + zion(iprot)*mp + zion(iprot)*me c 1 - bion(iprot)*mev2gr mion(iprot) = nion(iprot)*mn +zion(iprot)*mp -bion(iprot)*mev2gr wion(iprot) = avo * mion(iprot) wpart(iprot) = 2.0d0 ionam(iprot) = 'prot' c.. ionmax = ionmax + 1 ihe4 = ionmax aion(ihe4) = 4.0d0 nion(ihe4) = 2.0d0 zion(ihe4) = 2.0d0 bion(ihe4) = 28.29603d0 c mion(ihe4) = nion(ihe4)*mn + zion(ihe4)*mp + zion(ihe4)*me c 1 - bion(ihe4)*mev2gr mion(ihe4) = nion(ihe4)*mn + zion(ihe4)*mp - bion(ihe4)*mev2gr wion(ihe4) = avo * mion(ihe4) wpart(ihe4) = 1.0d0 ionam(ihe4) = 'he4' ionend = ionmax - 3 c..for up-arrow configurations else ionmax = i - 1 ionend = ionmax end if c..for either orientation, append energy, temperature, and denisty pointers iener = ionmax + 1 itemp = ionmax + 2 iden = ionmax + 3 ivelx = ionmax + 4 iposx = ionmax + 5 neqs = iposx ionam(iener) = 'ener ' ionam(itemp) = 'temp ' ionam(iden) = 'den ' ionam(ivelx) = 'velx ' ionam(iposx) = 'posx ' c..number of neutrons and do i = ionbeg,ionend nion(i) = aion(i) - zion(i) enddo c..mass of each isotope assuming fully ionized do i = ionbeg,ionend c mion(i) = nion(i)*mn + zion(i)*mp + zion(i)*me - bion(i)*mev2gr mion(i) = nion(i)*mn + zion(i)*mp - bion(i)*mev2gr enddo c..molar mass of each isotope do i = ionbeg,ionend wion(i) = avo * mion(i) enddo c..here is a common approximation do i=1,ionmax wion(i) = aion(i) enddo c..read data for electron capture on proton and positron capture on neutron c read(fil14,05) ((rrpen(j,i),i=1,7),j=1,6) c read(fil14,05) ((rrnep(j,i),i=1,7),j=1,6) read(fil14,05) ((xx,i=1,7),j=1,6) read(fil14,05) ((xx,i=1,7),j=1,6) c..build the links between the isotopes in the network c..before reading the weak reaction rates call naray c..read data for fuller weak rates. data is ordered in sequence of c..decreasing q-value for electron capture (in electron rest masses). icode c..keeps the matrix location of the isotope that is beta-decaying. data is c..tabular with 6 values of density and 7 of temperature. five quantities are c..tabulated: positron decay rate, effective electron capture ft value, beta c..decay rate, neutrino loss rate, and anti-neutrino loss rate. 140 read(fil14,06) nz,nn,qful if (nz.eq.99) go to 190 c..see if the isotope is in the network do i=ionbeg,ionend if (int(zion(i)) .eq. nz .and. int(nion(i)) .eq. nn 1 .and. nrr(2,i) .ne. 0) then nful = nful + 1 if (qful .gt. -1.0) nfulnot = nfulnot + 1 read(fil14,08) ((datful(nful,j,k),j=1,6),k=1,7) read(fil14,08) ((datful(nfuldim+nful,j,k),j=1,6),k=1,7) read(fil14,08) ((datful(2*nfuldim+nful,j,k),j=1,6),k=1,7) read(fil14,08) ((datful(3*nfuldim+nful,j,k),j=1,6),k=1,7) read(fil14,08) ((datful(4*nfuldim+nful,j,k),j=1,6),k=1,7) icode(nful) = i icode2(i) = nful qn(nful) = qful goto 140 end if enddo c..didn't find the isotope, or isotope in list but no link, still do the read do mm=1,35 read(fil14,07) xx enddo goto 140 c..finally all done reading the nuclear reaction rate file 190 continue close(unit=fil14) c..set the isotope names and pointers c write(6,*) ' ' c write(6,*) ' using isotopes:' c write(6,*) ' i z a name' if (.not.downarrow) then i = ineut inta = aion(i) intz = zion(i) c write(6,09) ineut,intz,inta,ineut,ionam(i) i = iprot inta = aion(i) intz = zion(i) c write(6,09) iprot,intz,inta,ionam(i) i = ihe4 inta = aion(i) intz = zion(i) c write(6,09) ihe4,intz,inta,ionam(i) endif do i=ionbeg,ionend inta = aion(i) intz = zion(i) if (intz .ge. 1 .and. intz .le. zmax) then if (inta .ge. aidmin(intz) .and. inta .le. aidmax(intz)) then do ii = aidmin(intz),aidmax(intz) if (ii .eq. inta) then write(string,01) zsymb(intz),inta call sqeeze(string) ionam(i) = string c write(6,09) i,intz,inta,ionam(i) end if enddo else write(6,*) ' bad aion',inta,' in routine init_torch' write(6,*) ' zion=',intz write(6,*) ' amin=',aidmin(intz),' amax=',aidmax(intz) stop 'error: bad inta in routine init_torch' end if else write(6,*) 'bad zion',intz,' in routine init_torch' write(6,*) 'inta =',inta,' zmax=', zmax stop 'error: bad intz in routine init_torch' end if enddo if (downarrow) then i = ineut inta = aion(i) intz = zion(i) c write(6,09) ionmax-2,intz,inta,ionam(i) i = iprot inta = aion(i) intz = zion(i) c write(6,09) ionmax-1,intz,inta,ionam(i) i = ihe4 inta = aion(i) intz = zion(i) c write(6,09) ionmax,intz,inta,ionam(i) endif c..check some things c do i=1,ionmax c write(6,888) ionam(i), c 1 int(zion(i)),int(nion(i)),int(aion(i)), c 1 mion(i)*avo,(mion(i)*avo-aion(i)), c 2 (mion(i)*avo-aion(i))/(mev2gr*avo) c 888 format(1x,a,3i4,1p5e18.10) c enddo c..set the id numbers of certain key isotopes do i=ionbeg,ionend if (ionam(i) .eq. 'h2 ') then ih2 = i else if (ionam(i) .eq. 'h3 ') then ih3 = i else if (ionam(i) .eq. 'he3 ') then ihe3 = i else if (ionam(i) .eq. 'li6 ') then ili6 = i else if (ionam(i) .eq. 'li7 ') then ili7 = i else if (ionam(i) .eq. 'li8 ') then ili8 = i else if (ionam(i) .eq. 'be7 ') then ibe7 = i else if (ionam(i) .eq. 'be9 ') then ibe9 = i else if (ionam(i) .eq. 'b8 ') then ib8 = i else if (ionam(i) .eq. 'b9 ') then ib9 = i else if (ionam(i) .eq. 'b10 ') then ib10 = i else if (ionam(i) .eq. 'b11 ') then ib11 = i else if (ionam(i) .eq. 'c11 ') then ic11 = i else if (ionam(i) .eq. 'c12 ') then ic12 = i else if (ionam(i) .eq. 'c13 ') then ic13 = i else if (ionam(i) .eq. 'c14 ') then ic14 = i else if (ionam(i) .eq. 'n12 ') then in12 = i else if (ionam(i) .eq. 'n13 ') then in13 = i else if (ionam(i) .eq. 'n14 ') then in14 = i else if (ionam(i) .eq. 'n15 ') then in15 = i else if (ionam(i) .eq. 'o14 ') then io14 = i else if (ionam(i) .eq. 'o15 ') then io15 = i else if (ionam(i) .eq. 'o16 ') then io16 = i else if (ionam(i) .eq. 'o17 ') then io17 = i else if (ionam(i) .eq. 'o18 ') then io18 = i else if (ionam(i) .eq. 'f17 ') then if17 = i else if (ionam(i) .eq. 'f18 ') then if18 = i else if (ionam(i) .eq. 'f19 ') then if19 = i else if (ionam(i) .eq. 'ne18 ') then ine18 = i else if (ionam(i) .eq. 'ne19 ') then ine19 = i else if (ionam(i) .eq. 'ne20 ') then ine20 = i else if (ionam(i) .eq. 'ne21 ') then ine21 = i else if (ionam(i) .eq. 'ne22 ') then ine22 = i else if (ionam(i) .eq. 'na20 ') then ina20 = i else if (ionam(i) .eq. 'na21 ') then ina21 = i else if (ionam(i) .eq. 'na22 ') then ina22 = i else if (ionam(i) .eq. 'na23 ') then ina23 = i else if (ionam(i) .eq. 'mg22 ') then img22 = i else if (ionam(i) .eq. 'mg23 ') then img23 = i else if (ionam(i) .eq. 'mg24 ') then img24 = i else if (ionam(i) .eq. 'mg25 ') then img25 = i else if (ionam(i) .eq. 'mg26 ') then img26 = i else if (ionam(i) .eq. 'al25 ') then ial25 = i else if (ionam(i) .eq. 'al26 ') then ial26 = i else if (ionam(i) .eq. 'al27 ') then ial27 = i else if (ionam(i) .eq. 'si27 ') then isi27 = i else if (ionam(i) .eq. 'si28 ') then isi28 = i else if (ionam(i) .eq. 'si29 ') then isi29 = i else if (ionam(i) .eq. 'si30 ') then isi30 = i else if (ionam(i) .eq. 'p30 ') then ip30 = i else if (ionam(i) .eq. 'p31 ') then ip31 = i else if (ionam(i) .eq. 's30 ') then is30 = i else if (ionam(i) .eq. 's31 ') then is31 = i else if (ionam(i) .eq. 's32 ') then is32 = i else if (ionam(i) .eq. 'cl35 ') then icl35 = i else if (ionam(i) .eq. 'ar36 ') then iar36 = i else if (ionam(i) .eq. 'k39 ') then ik39 = i else if (ionam(i) .eq. 'ca40 ') then ica40 = i else if (ionam(i) .eq. 'sc43 ') then isc43 = i else if (ionam(i) .eq. 'sc45 ') then isc45 = i else if (ionam(i) .eq. 'ti44 ') then iti44 = i else if (ionam(i) .eq. 'ti46 ') then iti46 = i else if (ionam(i) .eq. 'ti48 ') then iti48 = i else if (ionam(i) .eq. 'ti50 ') then iti50 = i else if (ionam(i) .eq. 'v46 ') then iv46 = i else if (ionam(i) .eq. 'v47 ') then iv47 = i else if (ionam(i) .eq. 'v48 ') then iv48 = i else if (ionam(i) .eq. 'v51 ') then iv51 = i else if (ionam(i) .eq. 'cr47 ') then icr47 = i else if (ionam(i) .eq. 'cr48 ') then icr48 = i else if (ionam(i) .eq. 'cr49 ') then icr49 = i else if (ionam(i) .eq. 'cr50 ') then icr50 = i else if (ionam(i) .eq. 'cr51 ') then icr51 = i else if (ionam(i) .eq. 'cr52 ') then icr52 = i else if (ionam(i) .eq. 'cr53 ') then icr53 = i else if (ionam(i) .eq. 'cr54 ') then icr54 = i else if (ionam(i) .eq. 'mn50 ') then imn50 = i else if (ionam(i) .eq. 'mn51 ') then imn51 = i else if (ionam(i) .eq. 'mn52 ') then imn52 = i else if (ionam(i) .eq. 'mn55 ') then imn55 = i else if (ionam(i) .eq. 'fe52 ') then ife52 = i else if (ionam(i) .eq. 'fe54 ') then ife54 = i else if (ionam(i) .eq. 'fe55 ') then ife55 = i else if (ionam(i) .eq. 'fe56 ') then ife56 = i else if (ionam(i) .eq. 'fe57 ') then ife57 = i else if (ionam(i) .eq. 'fe58 ') then ife58 = i else if (ionam(i) .eq. 'co54 ') then ico54 = i else if (ionam(i) .eq. 'co55 ') then ico55 = i else if (ionam(i) .eq. 'co56 ') then ico56 = i else if (ionam(i) .eq. 'co59 ') then ico59 = i else if (ionam(i) .eq. 'ni56 ') then ini56 = i else if (ionam(i) .eq. 'ni58 ') then ini58 = i else if (ionam(i) .eq. 'ni59 ') then ini59 = i else if (ionam(i) .eq. 'ni64 ') then ini64 = i else if (ionam(i) .eq. 'ni66 ') then ini66 = i else if (ionam(i) .eq. 'cu63 ') then icu63 = i else if (ionam(i) .eq. 'zn60 ') then izn60 = i else if (ionam(i) .eq. 'zn64 ') then izn64 = i end if enddo c..set reaction rate names and pointers nrat = 0 do j=ionbeg,ionend if (nrat+14 .gt. abigrat) stop 'nrat > abigrat in init_torch' c..set up the (n,g) and (g,n) names and a few reaction rate pointers k = nrr(1,j) if (k .gt. 0) then string = ionam(j)//'(n,g)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string string = ionam(k)//'(g,n)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ih2 .and. k .eq. ih3) then irdng = nrat-1 irtgn = nrat else if (j .eq. ili6 .and. k .eq. ili7) then irli6ng = nrat-1 irli7gn = nrat else if (j .eq. ili7 .and. k .eq. ili8) then irli7ng = nrat-1 irli8gn = nrat else if (j .eq. ine20 .and. k .eq. ine21) then irne20ng = nrat-1 irne21gn = nrat end if end if c..set up the (p,n) (n,p) beta- beta+ decay names and a few reaction rate pointers k = nrr(2,j) if (k .gt. 0) then string = ionam(j)//'(p,n)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string string = ionam(k)//'(n,p)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ih3 .and. k .eq. ihe3) then irtpn = nrat - 1 irhe3np = nrat else if (j .eq. ili7 .and. k .eq. ibe7) then irli7pn = nrat - 1 irbe7np = nrat else if (j .eq. ibe9 .and. k .eq. ib9) then irbe9pn = nrat - 1 irb9np = nrat else if (j .eq. ib11 .and. k .eq. ic11) then irb11pn = nrat - 1 irc11np = nrat else if (j .eq. ic13 .and. k .eq. in13) then irc13pn = nrat - 1 irn13np = nrat else if (j .eq. ic14 .and. k .eq. in14) then irc14pn = nrat - 1 irn13np = nrat else if (j .eq. in14 .and. k .eq. io14) then irn14pn = nrat - 1 iro14np = nrat else if (j .eq. in15 .and. k .eq. io15) then irn15pn = nrat - 1 iro15np = nrat else if (j .eq. if19 .and. k .eq. ine19) then irf19pn = nrat - 1 irne19np = nrat else if (j .eq. ine22 .and. k .eq. ina22) then irne22pn = nrat - 1 irna22np = nrat else if (j .eq. ina23 .and. k .eq. img23) then irna23pn = nrat - 1 irmg23np = nrat end if string = ionam(j)//'(n=>p)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ibe7 .and. k .eq. ili7) irbeec = nrat string = ionam(k)//'(p=>n)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string end if c..set up the (p,g) and (g,p) names and a few reaction rate pointers k = nrr(3,j) if (k .gt. 0) then string = ionam(j)//'(p,g)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string string = ionam(k)//'(g,p)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ih2 .and. k .eq. ihe3) then irdpg = nrat - 1 irhe3gp = nrat else if (j .eq. ili6 .and. k .eq. ibe7) then irli6pg = nrat - 1 irbe7gp = nrat else if (j .eq. ibe7 .and. k .eq. ib8) then irbe7pg = nrat - 1 irb8gp = nrat else if (j .eq. ibe9 .and. k .eq. ib10) then irbe9pg = nrat - 1 irb10gp = nrat else if (j .eq. ib10 .and. k .eq. ic11) then irb10pg = nrat - 1 irc11gp = nrat else if (j .eq. ib11 .and. k .eq. ic12) then irb11pg = nrat - 1 irc12gp = nrat else if (j .eq. ic11 .and. k .eq. in12) then irc11pg = nrat - 1 irn12gp = nrat else if (j .eq. ic14 .and. k .eq. in15) then irc14pg = nrat - 1 irn15gp = nrat else if (j .eq. in13 .and. k .eq. io14) then irn13pg = nrat - 1 iro14gp = nrat else if (j .eq. in14 .and. k .eq. io15) then irn14pg = nrat - 1 iro15gp = nrat else if (j .eq. in15 .and. k .eq. io16) then irn15pg = nrat - 1 iro16gp = nrat else if (j .eq. io16 .and. k .eq. if17) then iro16pg = nrat - 1 irf17gp = nrat else if (j .eq. io17 .and. k .eq. if18) then iro17pg = nrat - 1 irf18gp = nrat else if (j .eq. io18 .and. k .eq. if19) then iro18pg = nrat - 1 irf19gp = nrat else if (j .eq. if17 .and. k .eq. ine18) then irf17pg = nrat - 1 irne18gp = nrat else if (j .eq. if18 .and. k .eq. ine19) then irf18pg = nrat - 1 irne19gp = nrat else if (j .eq. if19 .and. k .eq. ine20) then irf19pg = nrat - 1 irne20gp = nrat else if (j .eq. ine19 .and. k .eq. ina20) then irne19pg = nrat - 1 irna20gp = nrat else if (j .eq. ine20 .and. k .eq. ina21) then irne20pg = nrat - 1 irna21gp = nrat else if (j .eq. ine21 .and. k .eq. ina22) then irne21pg = nrat - 1 irna22gp = nrat else if (j .eq. ine22 .and. k .eq. ina23) then irne22pg = nrat - 1 irna23gp = nrat else if (j .eq. ina21 .and. k .eq. img22) then irna21pg = nrat - 1 irmg22gp = nrat else if (j .eq. ina22 .and. k .eq. img23) then irna22pg = nrat - 1 irmg23gp = nrat else if (j .eq. ina23 .and. k .eq. img24) then irna23pg = nrat - 1 irmg24gp = nrat else if (j .eq. img24 .and. k .eq. ial25) then irmg24pg = nrat - 1 iral25gp = nrat else if (j .eq. img25 .and. k .eq. ial26) then irmg25pg = nrat - 1 iral26gp = nrat else if (j .eq. img26 .and. k .eq. ial27) then irmg26pg = nrat - 1 iral27gp = nrat else if (j .eq. ial25 .and. k .eq. isi26) then iral25pg = nrat - 1 irsi26gp = nrat else if (j .eq. ial26 .and. k .eq. isi27) then iral26pg = nrat - 1 irsi27gp = nrat else if (j .eq. ial27 .and. k .eq. isi28) then iral27pg = nrat - 1 irsi28gp = nrat else if (j .eq. isi27 .and. k .eq. ip28) then irsi27pg = nrat - 1 irp28gp = nrat else if (j .eq. isi28 .and. k .eq. ip29) then irsi28pg = nrat - 1 irp29gp = nrat else if (j .eq. isi29 .and. k .eq. ip30) then irsi29pg = nrat - 1 irp30gp = nrat else if (j .eq. isi30 .and. k .eq. ip31) then irsi30pg = nrat - 1 irp31gp = nrat end if end if c..set up the (a,p) and (p,a) names and a few reaction rate pointers k = nrr(4,j) if (k .gt. 0) then string = ionam(j)//'(a,p)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string string = ionam(k)//'(p,a)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ihe3 .and. k .eq. ili6) then irhe3ap = nrat - 1 irli6pa = nrat else if (j .eq. ili6 .and. k .eq. ibe9) then irli6ap = nrat - 1 irbe9pa = nrat else if (j .eq. ibe7 .and. k .eq. ib10) then irbe7ap = nrat - 1 irb10pa = nrat else if (j .eq. ib8 .and. k .eq. ic11) then irb8ap = nrat - 1 irc11pa = nrat else if (j .eq. ib11 .and. k .eq. ic14) then irb11ap = nrat - 1 irc14pa = nrat else if (j .eq. ic11 .and. k .eq. in14) then irc11ap = nrat - 1 irn14pa = nrat else if (j .eq. in13 .and. k .eq. io16) then irn13ap = nrat - 1 iro16pa = nrat else if (j .eq. in14 .and. k .eq. io17) then irn14ap = nrat - 1 iro17pa = nrat else if (j .eq. in15 .and. k .eq. io18) then irn15ap = nrat - 1 iro18pa = nrat else if (j .eq. io14 .and. k .eq. if17) then iro14ap = nrat - 1 irf17pa = nrat else if (j .eq. io15 .and. k .eq. if18) then iro15ap = nrat - 1 irf18pa = nrat else if (j .eq. io16 .and. k .eq. if19) then iro16ap = nrat - 1 irf19pa = nrat else if (j .eq. if17 .and. k .eq. ine20) then irf17ap = nrat - 1 irne20pa = nrat else if (j .eq. if19 .and. k .eq. ine22) then irf19ap = nrat - 1 irne22pa = nrat else if (j .eq. ine20 .and. k .eq. ina23) then irne20ap = nrat - 1 irna23pa = nrat else if (j .eq. ina21 .and. k .eq. img24) then irna21ap = nrat - 1 irmg24pa = nrat else if (j .eq. img24 .and. k .eq. ial27) then irmg24ap = nrat - 1 iral27pa = nrat else if (j .eq. img25 .and. k .eq. ial28) then irmg25ap = nrat - 1 iral28pa = nrat end if end if c..set up the (a,n) and (n,a) names and a few reaction rate pointers k = nrr(5,j) if (k .gt. 0) then string = ionam(j)//'(a,n)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string string = ionam(k)//'(n,a)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ih3 .and. k .eq. ili6) then irtan = nrat irli6na = nrat else if (j .eq. ili7 .and. k .eq. ib10) then irli7an = nrat - 1 irb10na = nrat else if (j .eq. ibe9 .and. k .eq. ic12) then irbe9an = nrat - 1 irc12na = nrat else if (j .eq. ib10 .and. k .eq. in13) then irb10an = nrat - 1 irn13na = nrat else if (j .eq. ib11 .and. k .eq. in14) then irb11an = nrat - 1 irn14na = nrat else if (j .eq. ic12 .and. k .eq. io15) then irc12an = nrat - 1 iro15na = nrat else if (j .eq. ic13 .and. k .eq. io16) then irc13an = nrat - 1 iro16na = nrat else if (j .eq. in14 .and. k .eq. if17) then irn14an = nrat - 1 irf17na = nrat else if (j .eq. in15 .and. k .eq. if18) then irn15an = nrat - 1 irf18na = nrat else if (j .eq. io17 .and. k .eq. ine20) then iro17an = nrat - 1 irne20na = nrat else if (j .eq. io18 .and. k .eq. ine21) then iro18an = nrat - 1 irne21na = nrat else if (j .eq. if19 .and. k .eq. ina22) then irf19an = nrat - 1 irna22na = nrat else if (j .eq. ine21 .and. k .eq. img24) then irne21an = nrat - 1 irmg24na = nrat else if (j .eq. ine22 .and. k .eq. img25) then irne22an = nrat - 1 irmg25na = nrat else if (j .eq. img25 .and. k .eq. isi28) then irmg25an = nrat - 1 irsi28na = nrat else if (j .eq. img26 .and. k .eq. isi29) then irmg26an = nrat - 1 irsi29na = nrat else if (j .eq. ial27 .and. k .eq. ip30) then iral27an = nrat - 1 irp30na = nrat end if end if c..and the (a,g) and (g,a) names and a few reaction rate pointers k = nrr(6,j) if (k .gt. 0) then string = ionam(j)//'(a,g)'//ionam(k) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string string = ionam(k)//'(g,a)'//ionam(j) call sqeeze(string) nrat = nrat + 1 ratnam(nrat) = string if (j .eq. ih2 .and. k .eq. ili6) then irdag = nrat irli6ga = nrat else if (j .eq. ih3 .and. k .eq. ili7) then irtag = nrat - 1 irli7ga = nrat else if (j .eq. ihe3 .and. k .eq. ibe7) then irhe3ag = nrat - 1 irbe7ga = nrat else if (j .eq. ili6 .and. k .eq. ib10) then irli6ag = nrat - 1 irb10ga = nrat else if (j .eq. ili7 .and. k .eq. ib11) then irli7ag = nrat - 1 irb11ga = nrat else if (j .eq. ibe7 .and. k .eq. ic11) then irbe7ag = nrat - 1 irc11ga = nrat else if (j .eq. ic12 .and. k .eq. io16) then ircag = nrat - 1 iroga = nrat else if (j .eq. ic14 .and. k .eq. io18) then irc14ag = nrat - 1 iro18ga = nrat else if (j .eq. in14 .and. k .eq. if18) then irn14ag = nrat - 1 irf18ga = nrat else if (j .eq. in15 .and. k .eq. if19) then irn15ag = nrat - 1 irf19ga = nrat else if (j .eq. io14 .and. k .eq. ine18) then iro14ag = nrat - 1 irne18ga = nrat else if (j .eq. io15 .and. k .eq. ine19) then iro15ag = nrat - 1 irne19ga = nrat else if (j .eq. io16 .and. k .eq. ine20) then iroag = nrat - 1 irnega = nrat else if (j .eq. io17 .and. k .eq. ine21) then iro17ag = nrat - 1 irne21ga = nrat else if (j .eq. io18 .and. k .eq. ine22) then iro18ag = nrat - 1 irne22ga = nrat else if (j .eq. ine20 .and. k .eq. img24) then irneag = nrat - 1 irmgga = nrat else if (j .eq. ine21 .and. k .eq. img25) then irne21ag = nrat - 1 irmg25ga = nrat else if (j .eq. ine22 .and. k .eq. img26) then irne22ag = nrat - 1 irmg26ga = nrat else if (j .eq. img24 .and. k .eq. isi28) then irmgag = nrat - 1 irsiga = nrat else if (j .eq. img25 .and. k .eq. isi29) then irmg25ag = nrat - 1 irsi29ga = nrat else if (j .eq. img26 .and. k .eq. isi30) then irmg26ag = nrat - 1 irsi30ga = nrat else if (j .eq. isi28 .and. k .eq. is32) then irsiag = nrat - 1 irsga = nrat end if end if enddo c..make sure we have the space before adding in the hardcoded rates if (nrat + 67 .gt. abigrat) stop 'nrat > abigrat in init_torch' c..for p(e-,nu)n and n(e+,nub)p reactions nrat = nrat + 1 irpen = nrat ratnam(irpen) = 'rpen' nrat = nrat + 1 irnep = nrat ratnam(irnep) = 'rnep' c..c12 reactions, first half of triple alpha if (ic12 .ne. 0) then nrat = nrat + 1 ir3a = nrat ratnam(ir3a) = 'r3a' nrat = nrat + 1 irg3a = nrat ratnam(irg3a) = 'ral' c..c12+c12 reactions; must have ne20, na23, mg23 in the network if (ine20 .ne. 0 .and. ina23 .ne. 0 .and. img23 .ne. 0) then nrat = nrat + 1 ir1212n = nrat ratnam(ir1212n) = 'r1212n' nrat = nrat + 1 irmg23nc = nrat ratnam(irmg23nc) = 'rmg23nc' nrat = nrat + 1 ir1212p = nrat ratnam(ir1212p) = 'r1212p' nrat = nrat + 1 irna23pc = nrat ratnam(irna23pc) = 'rna23pc' nrat = nrat + 1 ir1212a = nrat ratnam(ir1212a) = 'r1212a' nrat = nrat + 1 irne20ac = nrat ratnam(irne20ac) = 'rne20ac' end if end if c..o16+o16 reactions; must have si28, p 30, p31 and s31 in the network if (io16 .ne. 0 .and. isi28 .ne. 0 .and. ip30 .ne. 0 .and. 1 ip31 .ne. 0 .and. is31 .ne. 0) then nrat = nrat + 1 ir1616n = nrat ratnam(nrat) = 'r1616n' nrat = nrat + 1 irs31no = nrat ratnam(irs31no) = 'rs31no' nrat = nrat + 1 ir1616p = nrat ratnam(ir1616p) = 'r1616p' nrat = nrat + 1 irp31po = nrat ratnam(irp31po) = 'rp31po' nrat = nrat + 1 ir1616a = nrat ratnam(ir1616a) = 'r1616a' nrat = nrat + 1 irsi28ao = nrat ratnam(irsi28ao) = 'rsi28ao' nrat = nrat + 1 ir1616d = nrat ratnam(ir1616d) = 'r1616d' nrat = nrat + 1 irp30do = nrat ratnam(irp30do) = 'rp30do' end if c..c12+o16 reactions; must have mg24, al27, si27 in the network if (ic12 .ne. 0 .and. io16 .ne. 0 .and. 1 img24 .ne. 0 .and. ial27 .ne. 0 .and. isi27 .ne. 0) then nrat = nrat + 1 ir1216n = nrat ratnam(ir1216n) = 'r1216n' nrat = nrat + 1 irsi27no = nrat ratnam(irsi27no) = 'rsi27no' nrat = nrat + 1 ir1216p = nrat ratnam(ir1216p) = 'r1216p' nrat = nrat + 1 iral27po = nrat ratnam(iral27po) = 'ral27po' nrat = nrat + 1 ir1216a = nrat ratnam(ir1216a) = 'r1216a' nrat = nrat + 1 irmg24ao = nrat ratnam(irmg24ao) = 'rmg24ao' end if c..if we have deuterium if (ih2 .ne. 0) then c..pp nrat = nrat + 1 irpp = nrat ratnam(irpp) = 'rpp' nrat = nrat + 1 irpep = nrat ratnam(irpep) = 'rpep' c..p(n,g)d nrat = nrat + 1 irpng = nrat ratnam(irpng) = 'rpng' nrat = nrat + 1 irdgn = nrat ratnam(irdgn) = 'rdgn' c..d(p,n)2p nrat = nrat + 1 irdpn = nrat ratnam(irdpn) = 'rdpn' nrat = nrat + 1 ir2pnp = nrat ratnam(ir2pnp) = 'r2pnp' c..d(d,g)he4 nrat = nrat + 1 irddg = nrat ratnam(irddg) = 'rddg' nrat = nrat + 1 irhe4gd = nrat ratnam(irhe4gd) = 'rhe4gd' end if c..if we have tritium if (ih3 .ne. 0) then c..d(d,p)t nrat = nrat + 1 irddp = nrat ratnam(irddp) = 'rddp' nrat = nrat + 1 irtpd = nrat ratnam(irtpd) = 'rtpd' c..t(p,g)he4 nrat = nrat + 1 irh3pg = nrat ratnam(irh3pg) = 'rh3pghe4' nrat = nrat + 1 irhe4gp = nrat ratnam(irhe4gp) = 'rhe4gph3' c..t(d,n)he4 reaction nrat = nrat + 1 irtdn = nrat ratnam(irtdn) = 'rtdn' nrat = nrat + 1 irhe4nd = nrat ratnam(irhe4nd) = 'rhe4nd' c..t(t,2n)he4 nrat = nrat + 1 irtt2n = nrat ratnam(irtt2n) = 'rtt2n' nrat = nrat + 1 irhe42nt = nrat ratnam(irhe42nt) = 'rhe42nt' end if c..if we have he3 if (ihe3 .ne. 0) then c..he3(he3,2p)he4 nrat = nrat + 1 ir33 = nrat ratnam(ir33) = 'r33' nrat = nrat + 1 ir33inv = nrat ratnam(ir33inv) = 'r33inv' c..he3(p,e+nu)he4 nrat = nrat + 1 irhep = nrat ratnam(irhep) = 'rhep' c..he3(n,g)he4 nrat = nrat + 1 irhe3ng = nrat ratnam(irhe3ng) = 'rhe3ng' nrat = nrat + 1 irhe4gn = nrat ratnam(irhe4gn) = 'rhe4gnhe3' c..he3(d,p)he4 nrat = nrat + 1 irhe3dp = nrat ratnam(irhe3dp) = 'rhe3dp' nrat = nrat + 1 irhe4pd = nrat ratnam(irhe4pd) = 'rhe4pd' c..d(d,n)he3 nrat = nrat + 1 irddn = nrat ratnam(irddn) = 'rddn' nrat = nrat + 1 irhe3nd = nrat ratnam(irhe3nd) = 'rhe3nd' c..he3(t,d)he4 nrat = nrat + 1 irhe3td = nrat ratnam(irhe3td) = 'rhe3td' nrat = nrat + 1 irhe4dt = nrat ratnam(irhe4dt) = 'rhe4dt' c..he3(t,np)he4 nrat = nrat + 1 irhe3tnp = nrat ratnam(irhe3tnp) = 'rhe3tnp' end if c..li7 reactions if (ili7 .ne. 0) then c..li7(t,2n)2a nrat = nrat + 1 irli7t2n = nrat ratnam(irli7t2n) = 'rli7t2n' c..li7(p,g)be8 and li7(p,a)he4 nrat = nrat + 1 irli7pag = nrat ratnam(irli7pag) = 'rli7pag' nrat = nrat + 1 ir2he4ga = nrat ratnam(ir2he4ga) = 'r2he4ga' c..li7(d,n)2a nrat = nrat + 1 irli7dn = nrat ratnam(irli7dn) = 'rli7dn' c..li7(he3,np)2a nrat = nrat + 1 irli7he3np = nrat ratnam(irli7he3np) = 'rli7he3np' end if if (ibe7 .ne. 0) then c..be7(d,p)2a nrat = nrat + 1 irbe7dp = nrat ratnam(irbe7dp) = 'rbe7dp' c..be7(t,np)2a nrat = nrat + 1 irbe7tnp = nrat ratnam(irbe7tnp) = 'rbe7tnp' c..be7(he3,2p)2a nrat = nrat + 1 irbe7he32p = nrat ratnam(irbe7he32p) = 'rbe7he32p' end if if (ibe9 .ne. 0) then c..a(an,g)be9 nrat = nrat + 1 iraan = nrat ratnam(iraan) = 'raan' nrat = nrat + 1 irgaan = nrat ratnam(irgaan) = 'rgaan' c..be9(p,d)be8 =>2a nrat = nrat + 1 irbe9pd = nrat ratnam(irbe9pd) = 'rbe9pd' end if if (ib8 .ne. 0) then c..b8(p=>n)be8 =>2a reactions nrat = nrat + 1 irb8ep = nrat ratnam(irb8ep) = 'rb8ep' end if if (ib11 .ne. 0) then c..b11(p,a)be8 => 2a reactions nrat = nrat + 1 irb11pa = nrat ratnam(irb11pa) = 'rb11pa' nrat = nrat + 1 ir3ap = nrat ratnam(ir3ap) = 'r3ap' end if if (ic11 .ne. 0) then c..c11(n,a)be8 => 2a nrat = nrat + 1 irc11na = nrat ratnam(irc11na) = 'rc11na' end if c..say how many isotopes and rates are in this network c write(6,*) c write(6,*) 'ionmax=',ionmax,' nrates=',nrat c write(6,*) 'minimum size of cx array ',kmax c write(6,*) return end subroutine naray include 'implno.dek' include 'network.dek' c..this routine builds the nrr(7,i) and nrrneut(7,i) arrays, which specify c..the location of isotopes coupled to i by various reactions. c..the first index on nrr refers to reactions of the form c.. 1=ng 2=pn 3=pg 4=ap 5=an 6=ag 7=b- c.. c..while the first index on nrrneut refers to reactions of the form c.. 1=nu,e-,n 2=nu,e- 3=nu,e-,p 4=nu e+,n 5=nu,e+ 6=nu,e+,p c..declare integer i,k,n,kz,kn,jz(7),jn(7) c..initialize for nrr jz(1) = 0 jz(2) = 1 jz(3) = 1 jz(4) = 1 jz(5) = 2 jz(6) = 2 jz(7) = -1 jn(1) = 1 jn(2) = -1 jn(3) = 0 jn(4) = 2 jn(5) = 1 jn(6) = 2 jn(7) = 1 c..build nrr do i=ionbeg,ionend do n=1,7 nrr(n,i) = 0 kz = int(zion(i)) + jz(n) kn = int(nion(i)) + jn(n) do k=ionbeg,ionend if (kz.eq.int(zion(k)) .and. kn.eq.int(nion(k))) nrr(n,i)=k enddo enddo enddo c..initialize for nrrneut jz(1) = 1 jz(2) = 1 jz(3) = 0 jz(4) = -1 jz(5) = -1 jz(6) = -2 jz(7) = -2 jn(1) = -2 jn(2) = -1 jn(3) = -1 jn(4) = 0 jn(5) = 1 jn(6) = 1 jn(7) = -2 c..build nrrneut do i=ionbeg,ionend do n=1,7 nrrneut(n,i) = 0 kz = int(zion(i)) + jz(n) kn = int(nion(i)) + jn(n) do k=ionbeg,ionend if (kz.eq.int(zion(k)) .and. kn.eq.int(nion(k))) nrrneut(n,i)=k enddo enddo enddo return end subroutine azbar(xmass,aion,zion,wion,ionmax, 1 ymass,abar,zbar,wbar,ye,nxcess) include 'implno.dek' c..this routine calculates composition variables c..input: c..mass fractions = xmass(1:ionmax) dimensionless c..number of nucleons = aion(1:ionmax) dimensionless c..charge of nucleus = zion(1:ionmax) dimensionless c..atomic weight or molar mass = wion(1:ionmax) g/mole c..number of isotopes = ionmax c.. c..output: c..molar abundances = ymass(1:ionmax) mole/g c..mean number of nucleons = abar dimensionless c..mean nucleon charge = zbar dimensionless c..mean weight = wbar g/mole c..electron fraction = ye mole/g c..neutron excess = xcess c..declare the pass integer ionmax double precision xmass(ionmax),aion(ionmax),zion(ionmax), 1 wion(ionmax),ymass(ionmax),abar,zbar,wbar, 2 ye,nxcess c..local variables integer i double precision sum,sum1 c..molar abundances do i=1,ionmax ymass(i) = xmass(i)/wion(i) enddo c..mean molar mass sum = 0.0d0 do i=1,ionmax sum = sum + ymass(i) enddo wbar = 1.0d0/sum c..mean number of nucleons sum1 = 0.0d0 do i=1,ionmax sum1 = sum1 + aion(i)*ymass(i) enddo abar = wbar * sum1 c..mean charge sum = 0.0d0 do i=1,ionmax sum = sum + zion(i)*ymass(i) enddo zbar = wbar * sum c..electron fraction ye = sum c..neutron excess nxcess = sum1 - 2.0d0 * ye return end subroutine nse(tt,dd,yye,newguess,ipartition,icoulomb, 1 xmass_out,xmun,xmup,iprint) include 'implno.dek' include 'const.dek' include 'network.dek' c..this routine puts a chosen reaction network into its nse distribution. c..input: c..tt = temperature c..dd = density c..ye = electron mol number c..newguess = 0 = use the old initial guess c.. = 1 = a new initial guess is made c..ipartition = 0 = use temperature independent partiction functions c.. = 1 = use temperature dependent partition functions c..icoulomb = 0 = do not use coulomb corrections c.. = 1 = use coulomb corrections c..iprint = print flag c..output: c..xmass_out = output nse mass fractions c..xmun = chemical potential of neutrons c..xmup = chemical potential of protons c..declare the pass integer newguess,ipartition,icoulomb,iprint double precision tt,dd,yye,xmass_out(1),xmun,xmup c..communicate double precision temp,den,ye_want,beta,mu_c_p,yei(abignet), 1 mu_c(abignet) common /nsec1/ temp,den,ye_want,beta,mu_c_p,yei,mu_c c..locals external nsefunc logical check integer ntrial,nfev,ntaken,n,i parameter (ntrial = 200, n = 2) double precision x(n),amass,fac1,fac2,tolf,tolx,twopi, 1 dum,resid(n), 2 xne,ge,sqrtge,sqrtgi, 3 esqu,forthpi,third,fivth,a1,a2,a2inv, 4 rt3,half_rt3,a3 parameter (tolf = 1.0d-10, 1 tolx = 1.0d-12, 2 twopi = 2.0d0*pi, 3 esqu = qe*qe, 4 forthpi = 4.0d0 * pi/3.0d0, 5 third = 1.0d0/3.0d0, 6 fivth = 5.0d0/3.0d0, 7 a1 = -0.9052d0, 8 a2 = 0.6322d0, 9 a2inv = 1.0d0/a2, & rt3 = 1.7320508075688772d0, 1 half_rt3 = 0.5d0 * rt3) c..for the initial guess double precision hinv,mev2erg,mev2gr,a56,z56,n56,b56 parameter (hinv = 1.0d0/h, 1 mev2erg = ev2erg*1.0d6, 2 mev2gr = mev2erg/clight**2, 3 a56 = 56.0d0, 4 z56 = 28.0d0, 5 n56 = 28.0d0, 6 b56 = 4.8398d2) c..for the temperature dependent partition functions integer jd1,jd2,jd3,jd4,jd5,jxx double precision gi,g0,dg0,aa,t9,t9i,t92 integer ifirst data ifirst/1/ c..initialize a common quantity if (ifirst .eq. 1) then ifirst = 0 do i=1,ionmax yei(i) = zion(i)/aion(i) enddo end if c..fill the common block temp = tt den = dd ye_want = yye beta = 1.0d0/(kerg * temp) c..set the partition functions wpart(i). c..this assumes init_torch succesfully executed. c..wpart(i) = ground state patition function * g(i) c.. = as(jd1) * g(i) c..above t9=10 the fitting functions can go INF, so limit t9. t9 = min(temp * 1.0d-9, 10.0d0) t9i = 1.0d0/t9 t92 = t9*t9 do i=ionbeg,ionend jd1 = 5*(i-1) + 1 jd2 = jd1 + 1 jd3 = jd2 + 1 jd4 = jd3 + 1 jd5 = jd4 + 1 gi = 0.0d0 g0 = 1.0d0 if (as(jd2) .ne. 0.0) then aa = as(jd2)*t9i + as(jd3) + as(jd4)*t9 + as(jd5)*t92 gi = exp(aa) if (ist(i) .ne. 0) then do jxx=6*(i-1)+1,6*(i-1)+2*ist(i)-1,2 aa = gs(jxx+1) * exp(-gs(jxx)*t9i) g0 = g0 + aa enddo end if end if gi = g0 + gi c..using ground state spins only if (ipartition .eq. 0) then wpart(i) = as(jd1) c..or with temperature dependence else if (ipartition .eq. 1) then wpart(i) = as(jd1) * gi else stop 'unknown ipartition function value' end if enddo c..set the coulomb corrections mu_c_p = 0.0d0 do i=i,ionmax mu_c(i) = 0.0d0 enddo c..calder et al, apj 656 313 2007, eq a1 c..number density of free electrons from matter if (icoulomb .eq. 1) then xne = ye_want * avo * den ge = esqu * beta * (forthpi * xne)**third a3 = -0.5d0*sqrt(3.0d0) - a1 / sqrt(a2) do i=1,ionmax gi = zion(i)**(fivth) * ge sqrtgi = sqrt(gi) mu_c(i) = a1*(sqrt(gi*(a2+gi)) 1 - a2*log(sqrt(gi*a2inv) + sqrt(1.0d0 + gi*a2inv)) 2 + 2.0d0*a3*(sqrtgi - atan(sqrtgi))) enddo mu_c_p = mu_c(iprot) end if c write(6,*) mu_c_p,mu_c(iprot) c..here is an initial guess for the neutron and proton chemical potentials, c..(x1) and (x2) respectively. obtained by setting xmass(ini56) = 1, c..setting mup = mun, and inverting the saha equation. c..here we use pure ni56 as it emprically appears to be a c..robust guess for all temp, rho, ye combinations. if (newguess .eq. 1) then newguess = 0 amass = n56*mn + z56*mp - b56*mev2gr fac1 = a56/(avo * den) fac2 = twopi/(beta*h) * amass*hinv fac2 = fac2 * sqrt(fac2) x(1) = -(log(fac1*fac2)/beta + b56*ev2erg*1.0d6)/a56 x(2) = x(1) end if c..root find on mass and charge conservation for c..the chemical potentials of protons and neutrons call xnewt_nse(ntrial,x,n,tolx,tolf,ntaken,check,nfev,nsefunc) c..be sure we converged if (check .or. ntaken .eq. ntrial) then write(6,*) write(6,*) 'check convergence of root finder' write(6,*) end if c..some optional diagnostics if (iprint .eq. 1) then write(6,*) write(6,110) 'iterations taken =',ntaken write(6,110) 'function evals =',nfev write(6,111) 'roots =',x(1),x(2) call nsefunc(dum,x,resid) write(6,111) 'mass conservation residual =',resid(1) write(6,111) 'charge conservation residual =',resid(2) 110 format(1x,a,i4) 111 format(1x,a,1p2e14.6) end if c..fill the output array using the converged values call nsefunc(dum,x,resid) c..bound the output nse abundances do i=1,ionmax xmass_out(i) = min(1.0d0,max(xmass_nse(i),1.0d-30)) enddo xmun = x(1) xmup = x(2) return end subroutine nsefunc(x,y,f) include 'implno.dek' include 'const.dek' include 'network.dek' c..this routine returns the root find functions. c..input is the point x and y a vector of the unknowns. c..output is the vector of root find functions f, which should be the c..zero vector upon convergence. c..y(1) is input as the neutron chemical potential c..y(2) is input as the proton chemical potential c..declare the pass double precision x,y(*),f(*) c..locals integer i,indx(abignet),j,ifirst double precision ye,mu,fac1,fac2,fac3,sum,sum2, 1 deninv,ww,hinv,twopih,mev2erg,mev2gr parameter (hinv = 1.0d0/h, 1 twopih = 2.0d0 * pi/h, 2 mev2erg = ev2erg*1.0d6, 3 mev2gr = mev2erg/clight**2) c..communicate double precision temp,den,ye_want,beta,mu_c_p,yei(abignet), 1 mu_c(abignet) common /nsec1/ temp,den,ye_want,beta,mu_c_p,yei,mu_c c..chemical potential and mass fraction of each isotope c..hartmann et al, apj 297 837 1985, eq 2 c..calder et al, apj 656 313 2007, eq a1 deninv = 1.0d0/den ww = twopih/beta do i=1,ionmax mu = nion(i)*y(1) + zion(i)*y(2) c amass(i) = nion(i)*mn + zion(i)*mp - bion(i)*mev2gr fac1 = mion(i) * deninv * wpart(i) fac2 = ww * mion(i) * hinv fac2 = fac2*sqrt(fac2) fac3 = exp( beta * (mu + bion(i)*mev2erg) 1 - mu_c(i) + zion(i)*mu_c_p) xmass_nse(i) = fac1 * fac2 * fac3 enddo c..sum the mass fractions in ascending order to minimize roundoff call indexx(ionmax,xmass_nse,indx) sum = 0.0d0 do i=1,ionmax sum = sum + xmass_nse(indx(i)) c sum = sum + xmass_nse(i) enddo c..sum the mass fractions to form ye c..this formulation assumes mion(i) = aion(i)*amu ye = 0.0d0 do i=1,ionmax j = indx(i) ye = ye + yei(j) * xmass_nse(j) ye = ye + yei(i) * xmass_nse(i) enddo c..this formulation does not assume mion(i) = aion(i)*amu c..seitenzahl et al 2007 sum2 = 0.0d0 do i=1,ionmax sum2 = sum2 + amu/mion(i) * ((ye_want - 1.0d0)*zion(i) 1 + ye_want * nion(i)) * xmass_nse(i) enddo c..mass and charge conservation are the requirements f(1) = sum - 1.0d0 c f(2) = ye - ye_want f(2) = sum2 return end subroutine nsejac(x,y,f,dfdy,n,np) include 'implno.dek' include 'const.dek' include 'network.dek' c..this routine returns the functions and the jacobian to do the root find on c..input is x, and y(n) a vector of the unknowns. output is f(n) c..and its jacobian dfdy(np,np). c..y(1) is the neutron chemical potential c..y(2) is the proton chemical potential c..declare the pass integer n,np double precision x,y(n),f(n),dfdy(np,np) c..locals integer indx(abignet),i,j double precision mu,mubn,mubp, 1 fac1,fac2,fac3,fac4,fac5, 2 xmbn(abignet),xmbp(abignet),sum,sumbn,sumbp, 3 ye,yebn,yebp,deninv,ww,sum2,sum2bn,sum2bp, 4 hinv,twopih,mev2erg,mev2gr parameter (hinv = 1.0d0/h, 1 twopih = 2.0d0 * pi/h, 2 mev2erg = ev2erg*1.0d6, 3 mev2gr = mev2erg/clight**2) c..communicate double precision temp,den,ye_want,beta,mu_c_p,yei(abignet), 1 mu_c(abignet) common /nsec1/ temp,den,ye_want,beta,mu_c_p,yei,mu_c c..chemical potential and mass fraction of each isotope c..hartmann et al, apj 297 837 1985, eq 2 c..calder et al, apj 656 313 2007, eq a1 deninv = 1.0d0/den ww = twopih/beta c..loop over isotopes do i=1,ionmax mu = nion(i) * y(1) + zion(i) * y(2) mubn = nion(i) mubp = zion(i) c amass(i) = nion(i)*mn + zion(i)*mp - bion(i)*mev2gr fac1 = mion(i) * deninv * wpart(i) fac2 = ww * mion(i) * hinv fac2 = fac2 * sqrt(fac2) fac3 = exp( beta * (mu + bion(i) * ev2erg * 1.0d6) 1 - mu_c(i) + zion(i)*mu_c_p) fac4 = fac1 * fac2 * fac3 xmass_nse(i) = fac4 xmbn(i) = fac4 * beta * mubn xmbp(i) = fac4 * beta * mubp enddo c..sum the mass fractions in ascending order to minimize roundoff call indexx(ionmax,xmass_nse,indx) sum = 0.0d0 sumbn = 0.0d0 sumbp = 0.0d0 do i=1,ionmax j = indx(i) sum = sum + xmass_nse(j) sumbn = sumbn + xmbn(j) sumbp = sumbp + xmbp(j) c sum = sum + xmass_nse(i) c sumbn = sumbn + xmbn(i) c sumbp = sumbp + xmbp(i) enddo c..sum the mass fractions to form ye c..this formulation assumes mion(i) = aion(i)*amu ye = 0.0d0 yebn = 0.0d0 yebp = 0.0d0 do i=1,ionmax j = indx(i) fac5 = yei(j) ye = ye + fac5 * xmass_nse(j) yebn = yebn + fac5 * xmbn(j) yebp = yebp + fac5 * xmbp(j) c fac5 = yei(i) c ye = ye + fac5 * xmass_nse(i) c yebn = yebn + fac5 * xmbn(i) c yebp = yebp + fac5 * xmbp(i) enddo c..this formulation does not assume mion(i) = aion(i)*amu c..seitenzahl et al 2007 sum2 = 0.0d0 sum2bn = 0.0d0 sum2bp = 0.0d0 do i=1,ionmax fac5 = amu/mion(i) * ((ye_want - 1.0d0)*zion(i) 1 + ye_want * nion(i)) sum2 = sum2 + fac5 * xmass_nse(i) sum2bn = sum2bn + fac5 * xmbn(i) sum2bp = sum2bp + fac5 * xmbp(i) enddo c..mass and charge conservation are the requirements f(1) = sum - 1.0d0 c f(2) = ye - ye_want f(2) = sum2 c..jacobian dfdy(1,1) = sumbn dfdy(1,2) = sumbp c dfdy(2,1) = yebn c dfdy(2,2) = yebp dfdy(2,1) = sum2bn dfdy(2,2) = sum2bp return end subroutine xnewt_nse(ntrial,x,n,tolx,tolf,ntaken,check,nfev,func) include 'implno.dek' c..given an initial guess x(1:n) for the root of n equations, this routine c..finds the root by a globally convergent newtons method. the vector of c..functions to be zeroed, called fvec(1:n) in the routine below, is c..returned by the user supplied routine func. the output quantity check c..is false on nomal return and true if the routine has converged to a c..local minimum of the function xfminx_nse. if so, try restarting from a c..different initial guess. c..np is the maximum number of equations n c..ntrial is the maximum number of iterations to try c..ntaken is the number of iterations done c..tolf sets the convergence on function values c..tolmin sets the criterion for deciding wheather spurious convergence to c.. a false minimum of xfminx_nse has occured c..tolx is the convergence criteria on deltax c..stpmx is the scaled maximum step length allowed in the line searches c..nfev is the number of function evaluations c..declare the pass external func logical check integer ntrial,n,ntaken,nfev double precision x(n),tolx,tolf c..common block communicates values from routine xfminx_nse integer nn,np parameter (np = 4) double precision fvec(np) common /newtnse/ fvec,nn c..locals integer i,its,j,indx(np) double precision tolmin,stpmx,d,den,f,fold,stpmax,sum,temp,test, 1 fjac(np,np),g(np),p(np),xold(np),xfminx_nse,dum parameter (tolmin = 1.0d-12, 1 stpmx = 2.0d0) c..initialize if (n .gt. np) stop 'n > np in routine xnewt' nn = n f = xfminx_nse(x,func) nfev = 1 ntaken = 0 c.. test for the initial guess being a root, using a more stringent tolf test = 0.0d0 do i=1,n if (abs(fvec(i)) .gt. test) test = abs(fvec(i)) enddo if (test .lt. 0.01*tolf) then check = .false. return end if c..get stpmax for the line search sum = 0.0d0 do i=1,n sum = sum + x(i)*x(i) enddo stpmax = stpmx * max(sqrt(sum),dfloat(n)) c..start of iteration loop; get the jacobian do its = 1, ntrial ntaken = its c..second order accurate numerical jacobian c call jac_nse(dum,x,fjac,n,n,np,np,func) c nfev = nfev + 2*n + 1 c..analytic jacobian call nsejac(dum,x,fvec,fjac,n,np) nfev = nfev + 1 c..compute grad f for the line searches do i=1,n sum = 0.0d0 do j=1,n sum = sum + fjac(j,i)*fvec(j) enddo g(i) = sum enddo c..store x, and f and form right hand sides do i=1,n xold(i) = x(i) enddo fold = f do i=1,n p(i) = -fvec(i) enddo c..solve the linear systems c write(6,*) 'xnewt nse' c write(6,112) x(1),x(2) c write(6,112) fvec(1),fvec(2) c write(6,112) fjac(1,1),fjac(1,2),fjac(2,1),fjac(2,2) c 112 format(1x,1p2e14.6) call ludcmp_nse(fjac,n,np,indx,d) call lubksb_nse(fjac,n,np,indx,p) c..line search returns new x and f c..it also gets fvec at the new x when it calls xfminx_nse call lnsrch_nse(n,xold,fold,g,p,x,f,stpmax,check,nfev,func) c write(6,112) x(1),x(2) c write(6,112) f c write(6,*) c..test for convergence on function value test = 0.0d0 do i=1,n if (abs(fvec(i)) .gt. test) test = abs(fvec(i)) enddo if (test .lt. tolf) then check = .false. return end if c..check for zero gradiant, i.e. spurious convergence if (check) then test = 0.0d0 den = max(f, 0.5d0 * n) do i=1,n temp = abs(g(i)) * max(abs(x(i)),1.0d0)/den if (temp .gt. test) test = temp enddo if (test .lt. tolmin) then check = .true. else check = .false. end if return end if c..test for convergence on deltax test = 0.0d0 do i=1,n temp = (abs(x(i)-xold(i)))/max(abs(x(i)),1.0d0) if (temp .gt. test) test = temp enddo if (test .lt. tolx) return c write(6,*) its,test c..back for another iteration enddo check = .true. return end subroutine lnsrch_nse(n,xold,fold,g,p,x,f,stpmax,check,nfev,func) include 'implno.dek' c..given an n dimensional point xold(1:n), the value of the function fold c..and the gradient g(1:n) at the point, and a direction p(1:n), this routine c..finds a new point x(1:n) along the direction of p from xold where the c..function xfminx_nse has decreased "sufficiently". the new function value is c..returned in f. stpmax is an input quanity that limits the length of the c..steps so that the function is not evaluated in regions where it is c..undefined or subject to overflow. p is usually the newton direction. the c..output quantity check is false on normal exit, and true when x is too c..close to xold. in a minimization routine, this usually signals c..convergence and can be ignored. however, in a root finding routine, the c..calling routine should check wheather the convergence is spurious. c..declare the pass external func logical check integer n,nfev double precision f,fold,stpmax,g(n),p(n),x(n),xold(n) c..locals integer i double precision xfminx_nse,a,alam,alam2,alamin,b,disc,f2,rhs1, 1 rhs2,slope,sum,temp,test,tmplam, 2 alf,tolx,alam_start parameter (alf = 1.0d-4, 1 tolx = 3.0d-12, 2 alam_start = 1.0d0) c..alf ensures sufficient decrease in the function value, tolx is the c..convergence criterion on deltax c..initialize and scale if the attempted step is too big check = .false. sum = 0.0d0 do i=1,n sum = sum + p(i)*p(i) enddo sum = sqrt(sum) if (sum .gt. stpmax) then do i=1,n p(i) = p(i) * stpmax/sum enddo end if slope = 0.0d0 do i=1,n slope = slope + g(i)*p(i) enddo if (slope .ge. 0.0) stop 'roundoff problem in lnsrch_nse' c..compute lambda_min test = 0.0d0 do i=1,n temp = abs(p(i))/max(abs(xold(i)),1.0d0) if (temp .gt. test) test = temp enddo alamin = tolx/test c..always try a full newton step, start of iteration loop alam = alam_start 1 continue do i=1,n x(i) = xold(i) + alam*p(i) c..for the nse problem make sure the neutron and proton c..chemical potentials are less than or equal to zero c..hmm for low ye (0.2) and high density (1e14), the neutron c..chemical potential does indeed go positive. let's try c..letting it go wherever it wants c if (x(i) .gt. 0.0) x(i) = xold(i) + 1.0d-4*alam*p(i) if (x(i) .gt. 0.0) x(i) = 1.0d-10 enddo f = xfminx_nse(x,func) nfev = nfev + 1 c..convergence on deltax, for root finding, the calling routine c..should verify the convergence if (alam .lt. alamin) then do i=1,n x(i) = xold(i) enddo check = .true. return c..sufficient function decrease else if (f .le. fold + alf*alam*slope) then return c..backtrack else if (alam .eq. alam_start) then tmplam = -slope / (2.0d0 * (f-fold-slope)) else rhs1 = f - fold - alam*slope rhs2 = f2 - fold - alam2*slope a = (rhs1/alam**2 - rhs2/alam2**2)/(alam-alam2) b = (-alam2*rhs1/alam**2 + alam*rhs2/alam2**2) / (alam-alam2) if (a .eq. 0.0) then tmplam = -slope/(2.0d0 * b) else disc = b*b - 3.0d0 * a * slope if (disc .lt. 0.0) then tmplam = 0.5d0 * alam else if (b .le. 0.0) then tmplam = (-b + sqrt(disc)) / (3.0d0 * a) else tmplam = -slope/(b + sqrt(disc)) end if end if if (tmplam .gt. 0.5d0*alam) tmplam = 0.5d0*alam end if end if c..store for the next trip through alam2 = alam f2 = f alam = max(tmplam, 0.1d0*alam) goto 1 end double precision function xfminx_nse(x,func) include 'implno.dek' c..returns f = 0.5 f dot f at x. func is a user supplied routine of the c..functions to be root found. c..declare the pass external func double precision x(*) c..locals integer i double precision sum,dum c..common block communicates values back to routine xnewt integer nn,np parameter (np = 4) double precision fvec(np) common /newtnse/ fvec,nn call func(dum,x,fvec) sum = 0.0d0 do i=1,nn sum = sum + fvec(i)*fvec(i) enddo xfminx_nse = 0.5d0 * sum return end subroutine jac_nse(x,y,dfdy,mcol,nrow,mmax,nmax,derivs) include 'implno.dek' c..this routine computes a second order accurate jacobian matrix c..of the function contained in the routine derivs. c.. c..input is the point x and the the vector y(nrow) at which to compute the c..jacobian dfdy(mcol,nrow). c.. c..uses 2*nrow + 1 function evaluations c..declare the pass external derivs integer mcol,nrow,mmax,nmax double precision x,y(nmax),dfdy(mmax,nmax) c..locals integer i,j,imax parameter (imax = 4) double precision fminus(imax),fplus(imax),rel,ax,temp,h,hinv parameter (rel = 3.162278d-8, 1 ax = 1.0d-16) c..check if (nrow .gt. imax) stop 'nrow > imax in jacobian2' c..for each row, get the right stepsize do j=1,nrow temp = y(j) h = rel * max(abs(y(j)),ax) y(j) = temp + h h = y(j) - temp call derivs(x,y,fplus) y(j) = temp temp = y(j) y(j) = temp - h h = temp - y(j) call derivs(x,y,fminus) y(j) = temp c..compute the jth row of the jacobian hinv = 1.0d0/(2.0d0 * h) do i=1,mcol dfdy(i,j) = (fplus(i) - fminus(i)) * hinv enddo enddo c..restore the original state call derivs(x,y,fplus) return end subroutine helmeos include 'implno.dek' include 'const.dek' include 'vector_eos.dek' c..given a temperature temp [K], density den [g/cm**3], and a composition c..characterized by abar and zbar, this routine returns most of the other c..thermodynamic quantities. of prime interest is the pressure [erg/cm**3], c..specific thermal energy [erg/gr], the entropy [erg/g/K], along with c..their derivatives with respect to temperature, density, abar, and zbar. c..other quantites such the normalized chemical potential eta (plus its c..derivatives), number density of electrons and positron pair (along c..with their derivatives), adiabatic indices, specific heats, and c..relativistically correct sound speed are also returned. c.. c..this routine assumes planckian photons, an ideal gas of ions, c..and an electron-positron gas with an arbitrary degree of relativity c..and degeneracy. interpolation in a table of the helmholtz free energy c..is used to return the electron-positron thermodynamic quantities. c..all other derivatives are analytic. c.. c..references: cox & giuli chapter 24 ; timmes & swesty apj 1999 c..declare integer i,j double precision temp,den,abar,zbar,ytot1,ye, 1 x,y,zz,zzi,deni,tempi,xni,dxnidd,dxnida, 2 dpepdt,dpepdd,deepdt,deepdd,dsepdd,dsepdt, 3 dpraddd,dpraddt,deraddd,deraddt,dpiondd,dpiondt, 4 deiondd,deiondt,dsraddd,dsraddt,dsiondd,dsiondt, 5 dse,dpe,dsp,kt,ktinv,prad,erad,srad,pion,eion, 6 sion,xnem,pele,eele,sele,pres,ener,entr,dpresdd, 7 dpresdt,denerdd,denerdt,dentrdd,dentrdt,cv,cp, 8 gam1,gam2,gam3,chit,chid,nabad,sound,etaele, 9 detadt,detadd,xnefer,dxnedt,dxnedd,s double precision pgas,dpgasdd,dpgasdt,dpgasda,dpgasdz, 1 egas,degasdd,degasdt,degasda,degasdz, 2 sgas,dsgasdd,dsgasdt,dsgasda,dsgasdz, 3 cv_gas,cp_gas,gam1_gas,gam2_gas,gam3_gas, 4 chit_gas,chid_gas,nabad_gas,sound_gas double precision sioncon,forth,forpi,kergavo,ikavo,asoli3,light2 parameter (sioncon = (2.0d0 * pi * amu * kerg)/(h*h), 1 forth = 4.0d0/3.0d0, 2 forpi = 4.0d0 * pi, 3 kergavo = kerg * avo, 4 ikavo = 1.0d0/kergavo, 5 asoli3 = asol/3.0d0, 6 light2 = clight * clight) c..for the abar derivatives double precision dpradda,deradda,dsradda, 1 dpionda,deionda,dsionda, 2 dpepda,deepda,dsepda, 3 dpresda,denerda,dentrda, 4 detada,dxneda c..for the zbar derivatives double precision dpraddz,deraddz,dsraddz, 1 dpiondz,deiondz,dsiondz, 2 dpepdz,deepdz,dsepdz, 3 dpresdz,denerdz,dentrdz, 4 detadz,dxnedz c..for the tables, in general integer imax,jmax c..normal table c parameter (imax = 211, jmax = 71) c..big table c parameter (imax = 261, jmax = 101) c..bigger table parameter (imax = 271, jmax = 101) double precision d(imax),t(jmax) c..for the helmholtz free energy tables double precision f(imax,jmax),fd(imax,jmax), 1 ft(imax,jmax),fdd(imax,jmax),ftt(imax,jmax), 2 fdt(imax,jmax),fddt(imax,jmax),fdtt(imax,jmax), 3 fddtt(imax,jmax) c..for the pressure derivative with density ables double precision dpdf(imax,jmax),dpdfd(imax,jmax), 1 dpdft(imax,jmax),dpdfdt(imax,jmax) c..for chemical potential tables double precision ef(imax,jmax),efd(imax,jmax), 1 eft(imax,jmax),efdt(imax,jmax) c..for the number density tables double precision xf(imax,jmax),xfd(imax,jmax), 1 xft(imax,jmax),xfdt(imax,jmax) c..for the interpolations integer iat,jat double precision tlo,thi,tstp,tstpi,dlo,dhi,dstp,dstpi, 1 tsav,dsav,free,df_d,df_t,df_dd,df_tt,df_dt double precision dth,dt2,dti,dt2i,dd,dd2,ddi,dd2i,xt,xd,mxt,mxd, 1 si0t,si1t,si2t,si0mt,si1mt,si2mt, 2 si0d,si1d,si2d,si0md,si1md,si2md, 3 dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, 4 dsi0d,dsi1d,dsi2d,dsi0md,dsi1md,dsi2md, 5 ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt, 6 ddsi0d,ddsi1d,ddsi2d,ddsi0md,ddsi1md,ddsi2md, 7 z,psi0,dpsi0,ddpsi0,psi1,dpsi1,ddpsi1,psi2, 8 dpsi2,ddpsi2,din,h5,fi(36), 9 xpsi0,xdpsi0,xpsi1,xdpsi1,h3, 1 w0t,w1t,w2t,w0mt,w1mt,w2mt, 2 w0d,w1d,w2d,w0md,w1md,w2md c..for storing the differences double precision dt_sav(jmax),dt2_sav(jmax), 1 dti_sav(jmax),dt2i_sav(jmax), 2 dd_sav(imax),dd2_sav(imax), 3 ddi_sav(imax),dd2i_sav(imax) c..for the uniform background coulomb correction double precision dsdd,dsda,lami,inv_lami,lamida,lamidd, 1 plasg,plasgdd,plasgdt,plasgda,plasgdz, 3 ecoul,decouldd,decouldt,decoulda,decouldz, 4 pcoul,dpcouldd,dpcouldt,dpcoulda,dpcouldz, 5 scoul,dscouldd,dscouldt,dscoulda,dscouldz, 6 a1,b1,c1,d1,e1,a2,b2,c2,third,esqu parameter (a1 = -0.898004d0, 1 b1 = 0.96786d0, 2 c1 = 0.220703d0, 3 d1 = -0.86097d0, 4 e1 = 2.5269d0, 5 a2 = 0.29561d0, 6 b2 = 1.9885d0, 7 c2 = 0.288675d0, 8 third = 1.0d0/3.0d0, 9 esqu = qe * qe) c..for initialization integer ifirst data ifirst/0/ c..quintic hermite polynomial statement functions c..psi0 and its derivatives psi0(z) = z**3 * ( z * (-6.0d0*z + 15.0d0) -10.0d0) + 1.0d0 dpsi0(z) = z**2 * ( z * (-30.0d0*z + 60.0d0) - 30.0d0) ddpsi0(z) = z* ( z*( -120.0d0*z + 180.0d0) -60.0d0) c..psi1 and its derivatives psi1(z) = z* ( z**2 * ( z * (-3.0d0*z + 8.0d0) - 6.0d0) + 1.0d0) dpsi1(z) = z*z * ( z * (-15.0d0*z + 32.0d0) - 18.0d0) +1.0d0 ddpsi1(z) = z * (z * (-60.0d0*z + 96.0d0) -36.0d0) c..psi2 and its derivatives psi2(z) = 0.5d0*z*z*( z* ( z * (-z + 3.0d0) - 3.0d0) + 1.0d0) dpsi2(z) = 0.5d0*z*( z*(z*(-5.0d0*z + 12.0d0) - 9.0d0) + 2.0d0) ddpsi2(z) = 0.5d0*(z*( z * (-20.0d0*z + 36.0d0) - 18.0d0) + 2.0d0) c..biquintic hermite polynomial statement function h5(i,j,w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md)= 1 fi(1) *w0d*w0t + fi(2) *w0md*w0t 2 + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt 3 + fi(5) *w0d*w1t + fi(6) *w0md*w1t 4 + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt 5 + fi(9) *w0d*w2t + fi(10) *w0md*w2t 6 + fi(11) *w0d*w2mt + fi(12) *w0md*w2mt 7 + fi(13) *w1d*w0t + fi(14) *w1md*w0t 8 + fi(15) *w1d*w0mt + fi(16) *w1md*w0mt 9 + fi(17) *w2d*w0t + fi(18) *w2md*w0t & + fi(19) *w2d*w0mt + fi(20) *w2md*w0mt 1 + fi(21) *w1d*w1t + fi(22) *w1md*w1t 2 + fi(23) *w1d*w1mt + fi(24) *w1md*w1mt 3 + fi(25) *w2d*w1t + fi(26) *w2md*w1t 4 + fi(27) *w2d*w1mt + fi(28) *w2md*w1mt 5 + fi(29) *w1d*w2t + fi(30) *w1md*w2t 6 + fi(31) *w1d*w2mt + fi(32) *w1md*w2mt 7 + fi(33) *w2d*w2t + fi(34) *w2md*w2t 8 + fi(35) *w2d*w2mt + fi(36) *w2md*w2mt c..cubic hermite polynomial statement functions c..psi0 & derivatives xpsi0(z) = z * z * (2.0d0*z - 3.0d0) + 1.0 xdpsi0(z) = z * (6.0d0*z - 6.0d0) c..psi1 & derivatives xpsi1(z) = z * ( z * (z - 2.0d0) + 1.0d0) xdpsi1(z) = z * (3.0d0*z - 4.0d0) + 1.0d0 c..bicubic hermite polynomial statement function h3(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) = 1 fi(1) *w0d*w0t + fi(2) *w0md*w0t 2 + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt 3 + fi(5) *w0d*w1t + fi(6) *w0md*w1t 4 + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt 5 + fi(9) *w1d*w0t + fi(10) *w1md*w0t 6 + fi(11) *w1d*w0mt + fi(12) *w1md*w0mt 7 + fi(13) *w1d*w1t + fi(14) *w1md*w1t 8 + fi(15) *w1d*w1mt + fi(16) *w1md*w1mt c..popular format statements 01 format(1x,5(a,1pe11.3)) 02 format(1x,a,1p4e16.8) 03 format(1x,4(a,1pe11.3)) 04 format(1x,4(a,i4)) c..do this stuff once if (ifirst .eq. 0) then ifirst = 1 open(unit=19,file='helm_table.dat',status='old') c..read the normal helmholtz free energy table c tlo = 4.0d0 c thi = 11.0d0 c tstp = (thi - tlo)/float(jmax-1) c tstpi = 1.0d0/tstp c dlo = -10.0d0 c dhi = 11.0d0 c dstp = (dhi - dlo)/float(imax-1) c dstpi = 1.0d0/dstp c..for the bigger table c tlo = 3.0d0 c thi = 13.0d0 c tstp = (thi - tlo)/float(jmax-1) c tstpi = 1.0d0/tstp c dlo = -12.0d0 c dhi = 14.0d0 c dstp = (dhi - dlo)/float(imax-1) c dstpi = 1.0d0/dstp c..for the bigger table tlo = 3.0d0 thi = 13.0d0 tstp = (thi - tlo)/float(jmax-1) tstpi = 1.0d0/tstp dlo = -12.0d0 dhi = 15.0d0 dstp = (dhi - dlo)/float(imax-1) dstpi = 1.0d0/dstp do j=1,jmax tsav = tlo + (j-1)*tstp t(j) = 10.0d0**(tsav) do i=1,imax dsav = dlo + (i-1)*dstp d(i) = 10.0d0**(dsav) read(19,*) f(i,j),fd(i,j),ft(i,j),fdd(i,j),ftt(i,j),fdt(i,j), 1 fddt(i,j),fdtt(i,j),fddtt(i,j) enddo enddo c..read the pressure derivative with density table do j=1,jmax do i=1,imax read(19,*) dpdf(i,j),dpdfd(i,j),dpdft(i,j),dpdfdt(i,j) enddo enddo c..read the electron chemical potential table do j=1,jmax do i=1,imax read(19,*) ef(i,j),efd(i,j),eft(i,j),efdt(i,j) enddo enddo c..read the number density table do j=1,jmax do i=1,imax read(19,*) xf(i,j),xfd(i,j),xft(i,j),xfdt(i,j) enddo enddo c..construct the temperature and density deltas and their inverses do j=1,jmax-1 dth = t(j+1) - t(j) dt2 = dth * dth dti = 1.0d0/dth dt2i = 1.0d0/dt2 dt_sav(j) = dth dt2_sav(j) = dt2 dti_sav(j) = dti dt2i_sav(j) = dt2i end do do i=1,imax-1 dd = d(i+1) - d(i) dd2 = dd * dd ddi = 1.0d0/dd dd2i = 1.0d0/dd2 dd_sav(i) = dd dd2_sav(i) = dd2 ddi_sav(i) = ddi dd2i_sav(i) = dd2i enddo close(unit=19) c write(6,*) c write(6,*) 'finished reading eos table' c write(6,04) 'imax=',imax,' jmax=',jmax c write(6,03) 'temp(1) =',t(1),' temp(jmax) =',t(jmax) c write(6,03) 'ye*den(1) =',d(1),' ye*den(imax) =',d(imax) c write(6,*) end if c..start of pipeline loop, normal executaion starts here eosfail = .false. do j=jlo_eos,jhi_eos if (temp_row(j) .le. 0.0) stop 'temp less than 0 in helmeos' if (den_row(j) .le. 0.0) stop 'den less than 0 in helmeos' temp = temp_row(j) den = den_row(j) abar = abar_row(j) zbar = zbar_row(j) ytot1 = 1.0d0/abar ye = max(1.0d-16,ytot1 * zbar) c..initialize deni = 1.0d0/den tempi = 1.0d0/temp kt = kerg * temp ktinv = 1.0d0/kt c..radiation section: prad = asoli3 * temp * temp * temp * temp dpraddd = 0.0d0 dpraddt = 4.0d0 * prad*tempi dpradda = 0.0d0 dpraddz = 0.0d0 erad = 3.0d0 * prad*deni deraddd = -erad*deni deraddt = 3.0d0 * dpraddt*deni deradda = 0.0d0 deraddz = 0.0d0 srad = (prad*deni + erad)*tempi dsraddd = (dpraddd*deni - prad*deni*deni + deraddd)*tempi dsraddt = (dpraddt*deni + deraddt - srad)*tempi dsradda = 0.0d0 dsraddz = 0.0d0 c..ion section: xni = avo * ytot1 * den dxnidd = avo * ytot1 dxnida = -xni * ytot1 pion = xni * kt dpiondd = dxnidd * kt dpiondt = xni * kerg dpionda = dxnida * kt dpiondz = 0.0d0 eion = 1.5d0 * pion*deni deiondd = (1.5d0 * dpiondd - eion)*deni deiondt = 1.5d0 * dpiondt*deni deionda = 1.5d0 * dpionda*deni deiondz = 0.0d0 c..sackur-tetrode equation for the ion entropy of c..a single ideal gas characterized by abar x = abar*abar*sqrt(abar) * deni/avo s = sioncon * temp z = x * s * sqrt(s) y = log(z) sion = (pion*deni + eion)*tempi + kergavo * ytot1 * y dsiondd = (dpiondd*deni - pion*deni*deni + deiondd)*tempi 1 - kergavo * deni * ytot1 dsiondt = (dpiondt*deni + deiondt)*tempi - 1 (pion*deni + eion) * tempi*tempi 2 + 1.5d0 * kergavo * tempi*ytot1 x = avo*kerg/abar dsionda = (dpionda*deni + deionda)*tempi 1 + kergavo*ytot1*ytot1* (2.5d0 - y) dsiondz = 0.0d0 c..electron-positron section: c..assume complete ionization xnem = xni * zbar c..enter the table with ye*den din = ye*den c..bomb proof the input if (temp .gt. t(jmax)) then write(6,01) 'temp=',temp,' t(jmax)=',t(jmax) write(6,*) 'temp too hot, off grid' write(6,*) 'setting eosfail to true and returning' eosfail = .true. return end if if (temp .lt. t(1)) then write(6,01) 'temp=',temp,' t(1)=',t(1) write(6,*) 'temp too cold, off grid' write(6,*) 'setting eosfail to true and returning' eosfail = .true. return end if if (din .gt. d(imax)) then write(6,01) 'den*ye=',din,' d(imax)=',d(imax) write(6,*) 'ye*den too big, off grid' write(6,*) 'setting eosfail to true and returning' eosfail = .true. return end if if (din .lt. d(1)) then write(6,01) 'ye*den=',din,' d(1)=',d(1) write(6,*) 'ye*den too small, off grid' write(6,*) 'setting eosfail to true and returning' eosfail = .true. return end if c..hash locate this temperature and density jat = int((log10(temp) - tlo)*tstpi) + 1 jat = max(1,min(jat,jmax-1)) iat = int((log10(din) - dlo)*dstpi) + 1 iat = max(1,min(iat,imax-1)) c..access the table locations only once fi(1) = f(iat,jat) fi(2) = f(iat+1,jat) fi(3) = f(iat,jat+1) fi(4) = f(iat+1,jat+1) fi(5) = ft(iat,jat) fi(6) = ft(iat+1,jat) fi(7) = ft(iat,jat+1) fi(8) = ft(iat+1,jat+1) fi(9) = ftt(iat,jat) fi(10) = ftt(iat+1,jat) fi(11) = ftt(iat,jat+1) fi(12) = ftt(iat+1,jat+1) fi(13) = fd(iat,jat) fi(14) = fd(iat+1,jat) fi(15) = fd(iat,jat+1) fi(16) = fd(iat+1,jat+1) fi(17) = fdd(iat,jat) fi(18) = fdd(iat+1,jat) fi(19) = fdd(iat,jat+1) fi(20) = fdd(iat+1,jat+1) fi(21) = fdt(iat,jat) fi(22) = fdt(iat+1,jat) fi(23) = fdt(iat,jat+1) fi(24) = fdt(iat+1,jat+1) fi(25) = fddt(iat,jat) fi(26) = fddt(iat+1,jat) fi(27) = fddt(iat,jat+1) fi(28) = fddt(iat+1,jat+1) fi(29) = fdtt(iat,jat) fi(30) = fdtt(iat+1,jat) fi(31) = fdtt(iat,jat+1) fi(32) = fdtt(iat+1,jat+1) fi(33) = fddtt(iat,jat) fi(34) = fddtt(iat+1,jat) fi(35) = fddtt(iat,jat+1) fi(36) = fddtt(iat+1,jat+1) c..various differences xt = max( (temp - t(jat))*dti_sav(jat), 0.0d0) xd = max( (din - d(iat))*ddi_sav(iat), 0.0d0) mxt = 1.0d0 - xt mxd = 1.0d0 - xd c..the six density and six temperature basis functions si0t = psi0(xt) si1t = psi1(xt)*dt_sav(jat) si2t = psi2(xt)*dt2_sav(jat) si0mt = psi0(mxt) si1mt = -psi1(mxt)*dt_sav(jat) si2mt = psi2(mxt)*dt2_sav(jat) si0d = psi0(xd) si1d = psi1(xd)*dd_sav(iat) si2d = psi2(xd)*dd2_sav(iat) si0md = psi0(mxd) si1md = -psi1(mxd)*dd_sav(iat) si2md = psi2(mxd)*dd2_sav(iat) c..derivatives of the weight functions dsi0t = dpsi0(xt)*dti_sav(jat) dsi1t = dpsi1(xt) dsi2t = dpsi2(xt)*dt_sav(jat) dsi0mt = -dpsi0(mxt)*dti_sav(jat) dsi1mt = dpsi1(mxt) dsi2mt = -dpsi2(mxt)*dt_sav(jat) dsi0d = dpsi0(xd)*ddi_sav(iat) dsi1d = dpsi1(xd) dsi2d = dpsi2(xd)*dd_sav(iat) dsi0md = -dpsi0(mxd)*ddi_sav(iat) dsi1md = dpsi1(mxd) dsi2md = -dpsi2(mxd)*dd_sav(iat) c..second derivatives of the weight functions ddsi0t = ddpsi0(xt)*dt2i_sav(jat) ddsi1t = ddpsi1(xt)*dti_sav(jat) ddsi2t = ddpsi2(xt) ddsi0mt = ddpsi0(mxt)*dt2i_sav(jat) ddsi1mt = -ddpsi1(mxt)*dti_sav(jat) ddsi2mt = ddpsi2(mxt) c ddsi0d = ddpsi0(xd)*dd2i_sav(iat) c ddsi1d = ddpsi1(xd)*ddi_sav(iat) c ddsi2d = ddpsi2(xd) c ddsi0md = ddpsi0(mxd)*dd2i_sav(iat) c ddsi1md = -ddpsi1(mxd)*ddi_sav(iat) c ddsi2md = ddpsi2(mxd) c..the free energy free = h5(iat,jat, 1 si0t, si1t, si2t, si0mt, si1mt, si2mt, 2 si0d, si1d, si2d, si0md, si1md, si2md) c..derivative with respect to density df_d = h5(iat,jat, 1 si0t, si1t, si2t, si0mt, si1mt, si2mt, 2 dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md) c..derivative with respect to temperature df_t = h5(iat,jat, 1 dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, 2 si0d, si1d, si2d, si0md, si1md, si2md) c..derivative with respect to density**2 c df_dd = h5(iat,jat, c 1 si0t, si1t, si2t, si0mt, si1mt, si2mt, c 2 ddsi0d, ddsi1d, ddsi2d, ddsi0md, ddsi1md, ddsi2md) c..derivative with respect to temperature**2 df_tt = h5(iat,jat, 1 ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, 2 si0d, si1d, si2d, si0md, si1md, si2md) c..derivative with respect to temperature and density df_dt = h5(iat,jat, 1 dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, 2 dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md) c..now get the pressure derivative with density, chemical potential, and c..electron positron number densities c..get the interpolation weight functions si0t = xpsi0(xt) si1t = xpsi1(xt)*dt_sav(jat) si0mt = xpsi0(mxt) si1mt = -xpsi1(mxt)*dt_sav(jat) si0d = xpsi0(xd) si1d = xpsi1(xd)*dd_sav(iat) si0md = xpsi0(mxd) si1md = -xpsi1(mxd)*dd_sav(iat) c..derivatives of weight functions dsi0t = xdpsi0(xt)*dti_sav(jat) dsi1t = xdpsi1(xt) dsi0mt = -xdpsi0(mxt)*dti_sav(jat) dsi1mt = xdpsi1(mxt) dsi0d = xdpsi0(xd)*ddi_sav(iat) dsi1d = xdpsi1(xd) dsi0md = -xdpsi0(mxd)*ddi_sav(iat) dsi1md = xdpsi1(mxd) c..look in the pressure derivative only once fi(1) = dpdf(iat,jat) fi(2) = dpdf(iat+1,jat) fi(3) = dpdf(iat,jat+1) fi(4) = dpdf(iat+1,jat+1) fi(5) = dpdft(iat,jat) fi(6) = dpdft(iat+1,jat) fi(7) = dpdft(iat,jat+1) fi(8) = dpdft(iat+1,jat+1) fi(9) = dpdfd(iat,jat) fi(10) = dpdfd(iat+1,jat) fi(11) = dpdfd(iat,jat+1) fi(12) = dpdfd(iat+1,jat+1) fi(13) = dpdfdt(iat,jat) fi(14) = dpdfdt(iat+1,jat) fi(15) = dpdfdt(iat,jat+1) fi(16) = dpdfdt(iat+1,jat+1) c..pressure derivative with density dpepdd = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 si0d, si1d, si0md, si1md) dpepdd = max(ye * dpepdd,1.0d-30) c..look in the electron chemical potential table only once fi(1) = ef(iat,jat) fi(2) = ef(iat+1,jat) fi(3) = ef(iat,jat+1) fi(4) = ef(iat+1,jat+1) fi(5) = eft(iat,jat) fi(6) = eft(iat+1,jat) fi(7) = eft(iat,jat+1) fi(8) = eft(iat+1,jat+1) fi(9) = efd(iat,jat) fi(10) = efd(iat+1,jat) fi(11) = efd(iat,jat+1) fi(12) = efd(iat+1,jat+1) fi(13) = efdt(iat,jat) fi(14) = efdt(iat+1,jat) fi(15) = efdt(iat,jat+1) fi(16) = efdt(iat+1,jat+1) c..electron chemical potential etaele etaele = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 si0d, si1d, si0md, si1md) c..derivative with respect to density x = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 dsi0d, dsi1d, dsi0md, dsi1md) detadd = ye * x c..derivative with respect to temperature detadt = h3(iat,jat, 1 dsi0t, dsi1t, dsi0mt, dsi1mt, 2 si0d, si1d, si0md, si1md) c..derivative with respect to abar and zbar detada = -x * din * ytot1 detadz = x * den * ytot1 c..look in the number density table only once fi(1) = xf(iat,jat) fi(2) = xf(iat+1,jat) fi(3) = xf(iat,jat+1) fi(4) = xf(iat+1,jat+1) fi(5) = xft(iat,jat) fi(6) = xft(iat+1,jat) fi(7) = xft(iat,jat+1) fi(8) = xft(iat+1,jat+1) fi(9) = xfd(iat,jat) fi(10) = xfd(iat+1,jat) fi(11) = xfd(iat,jat+1) fi(12) = xfd(iat+1,jat+1) fi(13) = xfdt(iat,jat) fi(14) = xfdt(iat+1,jat) fi(15) = xfdt(iat,jat+1) fi(16) = xfdt(iat+1,jat+1) c..electron + positron number densities xnefer = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 si0d, si1d, si0md, si1md) c..derivative with respect to density x = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 dsi0d, dsi1d, dsi0md, dsi1md) x = max(x,1.0d-30) dxnedd = ye * x c..derivative with respect to temperature dxnedt = h3(iat,jat, 1 dsi0t, dsi1t, dsi0mt, dsi1mt, 2 si0d, si1d, si0md, si1md) c..derivative with respect to abar and zbar dxneda = -x * din * ytot1 dxnedz = x * den * ytot1 c..the desired electron-positron thermodynamic quantities c..dpepdd at high temperatures and low densities is below the c..floating point limit of the subtraction of two large terms. c..since dpresdd doesn't enter the maxwell relations at all, use the c..bicubic interpolation done above instead of the formally correct expression x = din * din pele = x * df_d dpepdt = x * df_dt c dpepdd = ye * (x * df_dd + 2.0d0 * din * df_d) s = dpepdd/ye - 2.0d0 * din * df_d dpepda = -ytot1 * (2.0d0 * pele + s * din) dpepdz = den*ytot1*(2.0d0 * din * df_d + s) x = ye * ye sele = -df_t * ye dsepdt = -df_tt * ye dsepdd = -df_dt * x dsepda = ytot1 * (ye * df_dt * din - sele) dsepdz = -ytot1 * (ye * df_dt * den + df_t) eele = ye*free + temp * sele deepdt = temp * dsepdt deepdd = x * df_d + temp * dsepdd deepda = -ye * ytot1 * (free + df_d * din) + temp * dsepda deepdz = ytot1* (free + ye * df_d * den) + temp * dsepdz c..coulomb section: c..uniform background corrections only c..from yakovlev & shalybkov 1989 c..lami is the average ion seperation c..plasg is the plasma coupling parameter z = forth * pi s = z * xni dsdd = z * dxnidd dsda = z * dxnida lami = 1.0d0/s**third inv_lami = 1.0d0/lami z = -third * lami lamidd = z * dsdd/s lamida = z * dsda/s plasg = zbar*zbar*esqu*ktinv*inv_lami z = -plasg * inv_lami plasgdd = z * lamidd plasgda = z * lamida plasgdt = -plasg*ktinv * kerg plasgdz = 2.0d0 * plasg/zbar c..yakovlev & shalybkov 1989 equations 82, 85, 86, 87 if (plasg .ge. 1.0) then x = plasg**(0.25d0) y = avo * ytot1 * kerg ecoul = y * temp * (a1*plasg + b1*x + c1/x + d1) pcoul = third * den * ecoul scoul = -y * (3.0d0*b1*x - 5.0d0*c1/x 1 + d1 * (log(plasg) - 1.0d0) - e1) y = avo*ytot1*kt*(a1 + 0.25d0/plasg*(b1*x - c1/x)) decouldd = y * plasgdd decouldt = y * plasgdt + ecoul/temp decoulda = y * plasgda - ecoul/abar decouldz = y * plasgdz y = third * den dpcouldd = third * ecoul + y*decouldd dpcouldt = y * decouldt dpcoulda = y * decoulda dpcouldz = y * decouldz y = -avo*kerg/(abar*plasg)*(0.75d0*b1*x+1.25d0*c1/x+d1) dscouldd = y * plasgdd dscouldt = y * plasgdt dscoulda = y * plasgda - scoul/abar dscouldz = y * plasgdz c..yakovlev & shalybkov 1989 equations 102, 103, 104 else if (plasg .lt. 1.0) then x = plasg*sqrt(plasg) y = plasg**b2 z = c2 * x - third * a2 * y pcoul = -pion * z ecoul = 3.0d0 * pcoul/den scoul = -avo/abar*kerg*(c2*x -a2*(b2-1.0d0)/b2*y) s = 1.5d0*c2*x/plasg - third*a2*b2*y/plasg dpcouldd = -dpiondd*z - pion*s*plasgdd dpcouldt = -dpiondt*z - pion*s*plasgdt dpcoulda = -dpionda*z - pion*s*plasgda dpcouldz = -dpiondz*z - pion*s*plasgdz s = 3.0d0/den decouldd = s * dpcouldd - ecoul/den decouldt = s * dpcouldt decoulda = s * dpcoulda decouldz = s * dpcouldz s = -avo*kerg/(abar*plasg)*(1.5d0*c2*x-a2*(b2-1.0d0)*y) dscouldd = s * plasgdd dscouldt = s * plasgdt dscoulda = s * plasgda - scoul/abar dscouldz = s * plasgdz end if c..bomb proof x = prad + pion + pele + pcoul y = erad + eion + eele + ecoul z = srad + sion + sele + scoul if (x .le. 0.0 .or. y .le. 0.0 .or. z .le. 0.0) then c write(6,*) c write(6,*) 'coulomb corrections are causing a negative pressure' c write(6,*) 'setting all coulomb corrections to zero' c write(6,*) pcoul = 0.0d0 dpcouldd = 0.0d0 dpcouldt = 0.0d0 dpcoulda = 0.0d0 dpcouldz = 0.0d0 ecoul = 0.0d0 decouldd = 0.0d0 decouldt = 0.0d0 decoulda = 0.0d0 decouldz = 0.0d0 scoul = 0.0d0 dscouldd = 0.0d0 dscouldt = 0.0d0 dscoulda = 0.0d0 dscouldz = 0.0d0 end if c pcoul = 0.0d0 c dpcouldd = 0.0d0 c dpcouldt = 0.0d0 c dpcoulda = 0.0d0 c dpcouldz = 0.0d0 c ecoul = 0.0d0 c decouldd = 0.0d0 c decouldt = 0.0d0 c decoulda = 0.0d0 c decouldz = 0.0d0 c scoul = 0.0d0 c dscouldd = 0.0d0 c dscouldt = 0.0d0 c dscoulda = 0.0d0 c dscouldz = 0.0d0 c..sum all the gas components pgas = pion + pele + pcoul egas = eion + eele + ecoul sgas = sion + sele + scoul dpgasdd = dpiondd + dpepdd + dpcouldd dpgasdt = dpiondt + dpepdt + dpcouldt dpgasda = dpionda + dpepda + dpcoulda dpgasdz = dpiondz + dpepdz + dpcouldz degasdd = deiondd + deepdd + decouldd degasdt = deiondt + deepdt + decouldt degasda = deionda + deepda + decoulda degasdz = deiondz + deepdz + decouldz dsgasdd = dsiondd + dsepdd + dscouldd dsgasdt = dsiondt + dsepdt + dscouldt dsgasda = dsionda + dsepda + dscoulda dsgasdz = dsiondz + dsepdz + dscouldz c..add in radiation to get the total pres = prad + pgas ener = erad + egas entr = srad + sgas dpresdd = dpraddd + dpgasdd dpresdt = dpraddt + dpgasdt dpresda = dpradda + dpgasda dpresdz = dpraddz + dpgasdz denerdd = deraddd + degasdd denerdt = deraddt + degasdt denerda = deradda + degasda denerdz = deraddz + degasdz dentrdd = dsraddd + dsgasdd dentrdt = dsraddt + dsgasdt dentrda = dsradda + dsgasda dentrdz = dsraddz + dsgasdz c..for the gas c..the temperature and density exponents (c&g 9.81 9.82) c..the specific heat at constant volume (c&g 9.92) c..th