*************** program hms11 c ----- calculates Hydrogen Main Sequence models, uses subroutines: c sch envel core step rhs state opact nburn solve prout c ----- kept as fortran files: c sch11 envel11 core11 step11 rhs1 state1 opact1 nburn1 solve11 prout1 c INPUT (from the keyboard): c X, Z = hydrogen and heavy element content (by mass fraction) c fmlz = log10 (stellar mass) (M Sun) c tel = log10 (effective temperature) (K) c fll = log10 (surface luminosity) (L Sun) c tcl = log10 (central temperature) (K) c rhcl = log10 (central density) (g/cm**3) c iprint = print control for the integrations, iprint<0 implies none c nm = number of stellar model to be calculated c dm = step in log10 (mass) between stellar models c NOTICE: tel, fll, tcl, rhcl are guessed values of the boundary parameters c to be improved by the subroutine sch c OUTPUT: c summary of results onto file "hms11.re1" and onto the screen c more detailes (when iprint>0) onto a file "hms11.re2" and onto the screen common/const/sunm,sunl,sunr,sunt,cgrrad,cdpm,sunmr3,sunml * ,pi,g,c,fmod c -- set the values of mathematical and physical constants in the common/const/ sunm=1.99e33 sunl=3.83e33 sunr=6.96e10 sunt=5800.0 pi=3.1415927 g=6.67e-8 c=2.998e10 fmod=alog(10.0) cgrrad=(sunl/sunm)/(16*pi*c*g) cdpm=g/4/pi*(sunm/sunr**2)**2 sunmr3=sunm/(4*pi*sunr**3) sunml=sunm/sunl c -- end of constants ------------------------------------------------------- write(*,103) 103 format(1x,/, *' program "hms11.for" calculates a series of hydrogen', *' main sequence models',/ *' with hydrogen content = X, and heavy element content = Z;',/ *' results are stored in files "hms11.re1" and "hms11.re2" '/) c Changed status='new' to status='unknown' in open statements --JJG open(1,file='hms11.re1',status='unknown') open(2,file='hms11.re2',status='unknown') 1 continue c -- enter parameters for the stellar models from the keyboard -------------- write(*,104) 104 format(' enter chemical composition: X, Z = ',$) read(*,*)x,z if(x.lt.0.0)stop write(2,100)x,z 100 format(1x,/,1x,' results of program "hms11.for", X, Z =', * 2f8.5) write(1,101) 101 format(1x,' X Z log M log Te log L log Tc' *,' log rhc log R iter') 2 continue write(*,105) 105 format(' log (Mass of the first star / solar mass) = ',$) read(*,*)fmlz write(*,106) 106 format(' log (effective temperature in degrees K) = ',$) read(*,*)tel write(*,107) 107 format(' log (luminosity / solar luminosity) = ',$) read(*,*)fll write(*,108) 108 format(' log (central temperature in degrees K) = ',$) read(*,*)tcl write(*,109) 109 format(' log (central density in g/cm**3) = ',$) read(*,*)rhcl write(*,110) 110 format(' output control (integer) = ',$) read(*,*)iprint write(*,111) 111 format(' number of stellar models = ',$) read(*,*)nm if(nm.le.0)go to 1 if(nm.eq.1)go to 4 write(*,112) 112 format(' step in log (Mass) between the models = ',$) read(*,*)dm 4 continue write(*,102) do 3 im=1,nm fml=fmlz+dm*(im-1) starm=10.0**fml call sch(x,z,starm,tel,fll,tcl,rhcl,rl,2,iprint,iter) write(*,101) write(*,102)x,z,fml,tel,fll,tcl,rhcl,rl,iter write(1,102)x,z,fml,tel,fll,tcl,rhcl,rl,iter 102 format(1x,f6.3,f10.7,6f8.4,i4) 3 continue write(*,102) go to 2 end ***************** ******************** subroutine sch(x,z,starm,tel,fll,tcl,rhcl,rl,iout,iprint,iter) c ----- calculates chemically homogeneous stellar model in a thermal c ----- equilibrium, i.e. a model on a zero age main sequence c ----- fitting results of integrations of envelope and core at the c ----- fitting mass, fmf=0.3*M c INPUT: c x = hydrogen content (by mass fraction) c z = heavy element content (by mass fraction) c starm = total stellar mass (M Sun) c tel = log10 (effective temperature) (K) (a guess) c fll = log10 (surface luminosity) (L Sun) (a guess) c tcl = log10 (central temperature) (K) (a guess) c rhcl = log10 (central density) (g/cm**3) (a guess) c iout = number of the output unit for storing results on the disk c iprint = output control c OUTPUT: c tel, fll, tcl, rhcl = accurate values obtained from numerical iterations c rl = log10 (surface radius) (R Sun) c iter = number of iterations needed to obtain those results c if iprint>0 then some results of integrations are stored on the disk c NOTICE: this subroutine uses subroutines: envel, core, solve c ---------- auxillary variables: c zn = nitrogen content for CNO burning (by mass fraction) c fmf = mass at the fitting point (M Sun) c difac = the maximum differences allowed at the fitting point c del = perturbations of the logarithms of boundary parameters c itmax = the maximum number of iterations c iter = the iteration counter c difmax = the maximum difference calculated at the fitting point c delfit(i) = the calculated differences at the fitting point for c the four variables = xe(i)/xc(i)-1 , i=1,2,3,4 c xe(i) = results of the envelope integrations at the fitting point c xc(i) = results of the core integrations at the fitting point c deriv(i,j) = partial derivatives of the differences at the fitting c point with respect to the boundary parameters c delce(i) = corrections added to the boundary parameters common/const/sunm,sunl,sunr,sunt,cgrrad,cdpm,sunmr3,sunml * ,pi,g,c,fmod dimension xe(10),xc(10),delfit(4),delce(4),deriv(4,5) * ,xe1(10),xe2(10),xc3(10),xc4(10) write(iout,100)starm,tel,fll,tcl,rhcl,iprint 100 format(1x,' Mass, log Te, log L, log Tc, log rhc, iter =',/, * f12.4,4f8.4,i5) c ----- set the constants for the fitting and the iterations ---- zn=z*0.3 fmf=starm*0.3 difac=0.0001 del=0.0001 itmax=15 iter=0 if(iprint.ge.(-1))write(*,104) 104 format(1x,' iter log Te log L log Tc log rhc', *' max del b.p. max non-fit') 2 continue c ----- a new iteration ---------------------------------------------------- iter=iter+1 if(iter.gt.itmax)go to 6 call envel(tel,fll,rl,xe,x,z,zn,fmf,starm,iout,iprint) call core(tcl,rhcl,xc,x,z,zn,fmf,iout,iprint) c ----- see if the envelope and the core integrations fit difmax=0.0 do 3 i=1,4 delfit(i)=xe(i)/xc(i)-1.0 dif=abs(delfit(i)) if(difmax.lt.dif)difmax=dif 3 continue if(difmax.lt.difac)go to 5 c ----- the differences at the fitting point are too large to be acceptable, c ----- calculate derivatives of differences at the fitting point with respect c ----- to the boundary parameters: tel, fll, tcl, rhcl call envel(tel+del,fll,rl,xe1,x,z,zn,fmf,starm,iout,iprint) call envel(tel,fll+del,rl,xe2,x,z,zn,fmf,starm,iout,iprint) call core(tcl+del,rhcl,xc3,x,z,zn,fmf,iout,iprint) call core(tcl,rhcl+del,xc4,x,z,zn,fmf,iout,iprint) do 4 i=1,4 deriv(i,1)=((xe1(i)/xc(i)-1.0)-delfit(i))/del deriv(i,2)=((xe2(i)/xc(i)-1.0)-delfit(i))/del deriv(i,3)=((xe(i)/xc3(i)-1.0)-delfit(i))/del deriv(i,4)=((xe(i)/xc4(i)-1.0)-delfit(i))/del deriv(i,5)=delfit(i) 4 continue c ----- calculate the corrections to the boundary parameters call solve(deriv,delce,facdel) c ----- add the corrections to the boundary parameters: tel=tel+delce(1) fll=fll+delce(2) tcl=tcl+delce(3) rhcl=rhcl+delce(4) if(iprint.lt.(-1))go to 2 write(*,103)iter,tel,fll,tcl,rhcl,facdel,difmax write(iout,103)iter,tel,fll,tcl,rhcl,facdel,difmax 103 format(1x,i5,4f8.4,3x,2f12.6) c ----- end of an iteration ------------------------------------------------ go to 2 5 continue c ---- iterations converged, stellar envelope and core fit each other if(iprint.ge.0)write(iout,102)tel,fll,tcl,rhcl,rl,iter 102 format(1x,' log Te, log L, log Tc, log rhc, log R, iter =' * /5f10.5,i5) return 6 continue write(*,*)' iterations have not converged ' write(iout,*)' iterations have not converged ' iter=-iter return end ******************* ******************* subroutine envel(tel,fll,rl,xx,x,z,zn,fmf,starm,iout,iprint) c -- integrates stellar structure equations from the surface to the fitting c -- point for a chemically homogeneous star c INPUT: c tel = log10 (effective temperature) (K) c fll = log10 (surface luminosity) (L Sun) c x = hydrogen content (by mass fraction) c z = heavy element content (by mass fraction) c zn = nitrogen content for CNO burning (by mass fraction) c fmf = mass at the fitting point (M Sun) c starm = total stellar mass (M Sun) c iout = the number of the output unit for subroutine prout c iprint = print control, one out of "iprint" steps to be printed with c subroutine prout, if iprint<0 then no print c OUTPUT: c rl = log10 (surface radius) (R Sun) c xx(i), i=1,2,3,4,5,6,7,8 at the fitting point c xx(1) = temperature (K) c xx(2) = density (g/cm**3) c xx(3) = radius (R Sun) c xx(4) = luminosity (L Sun) c xx(5) = mass (M Sun) c xx(6) = x = hydrogen content (by mass fraction) c xx(7) = z = heavy element content (by mass fraction) c xx(8) = zn = nitrogen content for CNO burning (by mass fraction) c NOTICE: subroutine envel1 uses: c subroutine step to make integration steps c subroutine prout to output results onto the terminal and unit "iout" common/const/sunm,sunl,sunr,sunt,cgrrad,cdpm,sunmr3,sunml * ,pi,g,c,fmod dimension xx(10) c --- calculate the conditions at the photosphere: temperature (K), c --- luminosity (L Sun), and radius (R Sun) tph=10.0**tel flph=10.0**fll rph=sqrt(flph)/(tph/sunt)**2 rl=alog10(rph) c --- set all the variables at the outer boundary, at a very small density, c --- i.e. at a very small optical depth, where temperature is lower than c --- the photospheric by a factor 2**0.25; assume the atmosphere is very c --- thin geometrically, so the radius is the same as at the photosphere. xx(1)=tph/2.0**0.25 xx(2)=1.0e-12 xx(3)=rph xx(4)=flph xx(5)=starm xx(6)=x xx(7)=z xx(8)=zn if(iprint.gt.0)call prout(0,xx,iout,-1) if(iprint.gt.0)call prout(0,xx,iout,1) kmax=1000 ip=iprint c ----- integrate the structure equations inward to the fitting point do 1 k=1,kmax call step(xx,fmf) ip=ip-1 if(ip.ne.0)go to 2 call prout(k,xx,iout,1) ip=iprint 2 continue if(abs(xx(5)/fmf-1.0).lt.0.00001)go to 3 1 continue write(*,100) write(*,iout) 100 format(1x,' too many integration steps in subroutine envel') call prout(k,xx,iout,1) stop 3 continue c ------ integrations completed ----- if(iprint.gt.0)call prout(-k,xx,iout,1) return end ***************** **************** subroutine core(tcl,rhcl,xx,x,z,zn,fmf,iout,iprint) c -- integrates stellar structure equations from the center to the fitting c -- point for a chemically homogeneous star c INPUT: c rhcl = log10 (central density) (g/cm**3) c tcl = log10 (central temperature) (K) c x = hydrogen content (by mass fraction) c z = heavy element content (by mass fraction) c zn = nitrogen content for CNO burning (by mass fraction) c fmf = mass at the fitting point (M Sun) c iout = the number of the output unit for subroutine prout c iprint = print control, one out of "iprint" steps to be printed with c subroutine prout, if iprint<0 then no printout c OUTPUT: c xx(i), i=1,2,3,4,5,6,7,8 at the fitting point c xx(1) = temperature (K) c xx(2) = density (g/cm**3) c xx(3) = radius (R Sun) c xx(4) = luminosity (L Sun) c xx(5) = mass (M Sun) c xx(6) = x = hydrogen content (by mass fraction) c xx(7) = z = heavy element content (by mass fraction) c xx(8) = zn = nitrogen content for CNO burning (by mass fraction) c NOTICE: subroutine core uses: c subroutine nburn to calculate luminosity at the inner boundary c subroutine step to make integration steps c subroutine prout to output results onto the terminal and unit "iout" c --------------- variables from common/const/ used in this subroutine: c sunmr3 = ( solar mass ) / ( 4*pi*(solar radius)**3 ) (c.g.s.) c sunml = ( solar mass ) / ( solar luminosity ) (c.g.s.) c --------------- auxillary variables: c tc = central temperature (K) c rhc = central density (g/cm**3) c rc = radius of the central sphere / solar radius c flc = luminosity generated in the central sphere / solar luminosity c fmc = mass of the central sphere / solar mass c epsx = energy generation rate in hydrogen burning (erg/g/sec) common/const/sunm,sunl,sunr,sunt,cgrrad,cdpm,sunmr3,sunml * ,pi,g,c,fmod dimension xx(10) c --- calculate variables at the surface of the inner sphere, i.e. inner c --- boundary conditions (this section may be a separate subroutine center) tc=10.0**tcl rhc=10.0**rhcl fmc=fmf*0.0001 call nburn(rhc,tc,x,zn,epsx,dxt) rc=(fmc/rhc*sunmr3*3)**(1.0/3.0) flc=sunml*fmc*epsx xx(1)=tc xx(2)=rhc xx(3)=rc xx(4)=flc xx(5)=fmc xx(6)=x xx(7)=z xx(8)=zn c ----- all variables at the starting point are ready if(iprint.gt.0)call prout(0,xx,iout,-1) if(iprint.gt.0)call prout(0,xx,iout,1) kmax=1000 ip=iprint c --- integrate the equations out to the fitting point: do 1 k=1,kmax call step(xx,fmf) ip=ip-1 if(ip.ne.0)go to 2 call prout(k,xx,iout,1) ip=iprint 2 continue if(abs(xx(5)/fmf-1.0).lt.0.00001)go to 3 1 continue write(*,100) write(iout,100) 100 format(1x,' too many integration steps in the subroutine core1 ') call prout(k,xx,iout,1) stop 3 continue c ----- integrations completed ---- if(iprint.gt.0)call prout(-k,xx,iout,1) return end ***************** **************** subroutine step(xx,fmf) c -- makes one integration step for a stellar core or envelope with a second c -- order Runge-Kutta method, the steps are carried up to the point xx(5) = fmf c INPUT: c xx(1) = t temperature (K) c xx(2) = rh density (g/cm**3) c xx(3) = r radius (R Sun) c xx(4) = flr luminosity (L Sun) c xx(5) = fmr mass (M Sun) (independent variable) c fmf mass at the fitting point (M Sun) c OUTPUT: c xx(i), i=1,2,3,4,5 at the end of integration step c NOTICE: subroutine step uses: c subroutine rhs to calculate right hand sides of differential equations c ------------- auxillary variables: c acc(i), i=1,2,3,4,5 = maximum step size in all variables c h = integration step size c h2 = half of the integration step c yy(i) = d xx(i) / d xx(nv) as calculated by subroutine rhs1 c xi(i) = values of xx(i) variables at the middle of the integration step dimension xx(10),yy(10),xi(10),acc(5) data acc/0.05,0.15,0.05,0.2,0.2/ call rhs(xx,yy) c ------ estimate the integration step = h ----- hi=abs(yy(5)/acc(5)/xx(5)) do 3 i=1,4 h=abs(yy(i)/acc(i)/xx(i)) if(hi.lt.h)hi=h 3 continue h=1.0/hi if(xx(5).gt.fmf)h=-h if(((xx(5)-fmf)*(xx(5)+h-fmf)).lt.0.0)h=fmf-xx(5) h2=h*0.5 c ----- make the first half of the integration step ------- do 1 i=1,5 xi(i)=xx(i)+h2*yy(i) 1 continue xi(6)=xx(6) xi(7)=xx(7) xi(8)=xx(8) call rhs(xi,yy) c ----- make the whole integration step ------- do 2 i=1,5 xx(i)=xx(i)+h*yy(i) 2 continue return end **************** ************* subroutine rhs(xx,yy) c -- calculates right hand sides of the stellar structure equations c INPUT: c xx(1) = t (K) temperature c xx(2) = rh (g/cm**3) density c xx(3) = r (R Sun) radius c xx(4) = flr (L Sun) luminosity c xx(5) = fmr (M Sun) mass (independent variable) c xx(6) = x hydrogen content (by mass fraction) c xx(7) = z heavy element content (by mass fraction) c xx(8) = zn nitrogen content for CNO burning (by mass) c if xx(8) < 0.0 then there is no nuclear burning - this is c to be used for example in envelope integrations c OUTPUT: c yy(i) = d xx(i) / d xx(1) i=1,2,3,4,5 c --------- NOTICE: subroutine rhs uses: c subroutine state to calculate thermodynamic functions c subroutine opact to calculate opacity (total: radiative and conductive) c subroutine nburn to calculate energy generation rate in hydrogen burning c --------------- variables from common/const/ used by this subroutine: c cgrrad = (sunl/sunm)/(16*pi*c*G) c cdpm = (G*sunm**2)/(4*pi*sunr**4) c sunmr3 = (sunm)/(4*pi*sunr**3) c sunml = (sunm)/(sunl) c sunm, sunl, sunr are the solar mass, luminosity and radius c --------------- auxillary variables: c p = total pressure (c.g.s) c pt = d ln p / d ln t at constant density c pr = d ln p / d ln rh at constant temperature c prad = radiation pressure (c.g.s.) c grad = d ln t / d ln p at constant entropy (adiabatic temperature c gradient) c grrad = d ln t / d ln p radiative temperature gradient c grt = d ln t / d ln p temperature gradient in the star c grrh = d ln rh / d ln p logarithmic density gradient in the star c dpm = d ln p / d ( Mr / M Sun ) pressure gradient common/const/sunm,sunl,sunr,sunt,cgrrad,cdpm,sunmr3,sunml * ,pi,g,c,fmod dimension xx(10),yy(10) t=xx(1) rh=xx(2) r=xx(3) flr=xx(4) fmr=xx(5) x=xx(6) z=xx(7) zn=xx(8) epsx=0.0 call state(rh,t,x,z,p,pt,pr,pgas,prad,grad,qt,qr) call opact(rh,t,x,z,fkap) if(zn.gt.0.0)call nburn(rh,t,x,zn,epsx,dxt) grrad=cgrrad*fkap*p/prad*flr/fmr grt=grad if(grt.gt.grrad)grt=grrad grrh=(1-grt*pt)/pr dpm=-cdpm*fmr/r**4 yy(1)=grt*dpm*t/p yy(2)=grrh*dpm*rh/p yy(3)=sunmr3/rh/r**2 yy(4)=sunml*epsx yy(5)=1.0 return end ***************** ************** subroutine state(rh,t,x,z,p,pt,pr,pgas,prad,grad,qt,qr) c -- calculates equation of state and thermodynamic quantities c INPUT: c rh = density (g/cm**2) c t = temperature (K) c x = hydrogen content (by mass fraction) c z = heavy element content (by mass fraction; helium content: y=1-x-z) c OUTPUT: c p = total pressure (radiation+ions+electrons; dyn/cm**2) c pt = d log p / d log t (at constant density) c pr = d log p / d log rh (at constant temperature) c prad = radiation pressure c pgas = gas pressure (ions+electrons) c grad = adiabatic temperature gradient: (d log t / d log p) at constant c entropy c qt = t*cv, where "cv" is specific heat at constant volume (per gram) c qr = qt * ( d log t / d log rh) at constant entropy c ---------------- NOTICE: full ionization is assumed, but electrons may be c partly degenerate, partly relativistic; c all formulae are written so as to avoid overflows; c ------ auxiliary variables: c y = helium content (by mass) c fmi = mean molecular weight of ions (oxygen stands for heavy elements) c fme = mean molecular weight per free electron (full ionization) c pion = ion pressure (perfect gas law) c pedr = pressure of relativistically degenerate electrons c pednr = pressure of non-relativistically degenerate electrons c ped = combined pressure of degenerate electrons c pend = pressure on non-degenerate electrons (perfect gas law) c pe = combined electron pressure (partly relativistic and degenerate) c f = a factor indicating how relativistic are the electrons c pet = d log pe / d log t (at constant density) c per = d log pe / d log rh (at constant temperature) data fkh,fk1,fk2,arad/0.8251e8, 9.91e12, 1.231e15, 2.523e-15/ y=1-x-z fmi=1/(x+y/4+z/16) fme=2/(1+x) prad=arad*t**4 pion=fkh/fmi*rh*t pedr=fk2*(rh/fme)**(4.0/3.0) pednr=fk1*(rh/fme)**(5.0/3.0) ped=pednr/sqrt(1+(pednr/pedr)**2) pend=fkh/fme*rh*t pe=pend*sqrt(1+(ped/pend)**2) pgas=pion+pe p=prad+pgas f=5.0/3.0*(ped/pednr)**2+4.0/3.0*(ped/pedr)**2 pet=(pend/pe)**2 per=pet+(1-pet)*f pt=(4*prad+pion+pet*pe)/p pr=(pion+per*pe)/p qt=(12*prad+1.5*pion+pet/(f-1)*pe)/rh qr=p/rh*pt grad=1/(qt/qr*pr+pt) return end ************** ************ subroutine nburn(rh,t,x,zn,epsx,dxt) c - calculates hydrogen burning rates for SCH and HEN programs, Oct. 22, 1984 c INPUT: c rh = density (g/cm**3) c t = temperature (K) c x = hydrogen content (by mass fraction) c zn = nitrogen content (by mass fraction) c OUTPUT: c epsx = energy generation rate due to p-p and CNO hydrogen burning c (erg/g/sec) c dxt = hydrogen depletion rate (g/g/sec) (negative number) c ------ auxillary variables: c epscno = energy generation rate in CNO cycle c epspp = energy generation rate in the pp reaction epsx=1.0e-30 if(t.lt.1.0e6)go to 1 c ------ if temperature is below 10**6 K the burning rate is fixed at 1.0e-30 c ------ to avoid underflows t613=(t/1.0e6)**(1.0/3.0) t623=t613*t613 epscno=x*zn*rh*exp(64.24-152.313/t613)/t623 epspp=x*x*rh*exp(14.54-33.81/t613)/t623 epsx=epscno+epspp 1 continue dxt=-1.667e-19*epsx return end ************* *********** subroutine opact(rh,t,x,z,fkap) c -- calculates opacity (radiative and conductive) c INPUT: c rh = density (g/cm**3) c t = temperature (K) c x = hydrogen content (by mass fraction) c z = heavy element content (by mass fraction; helium: y=1-x-z) c OUTPUT: c fkap = opacity (cm**2/g) c ----------- auxiliary variables: c y = helium conetent (by mass) c zav = average charge of ions for conductive opacity c fke = electron scattering opacity c fkk = "Kramers" (i.e. bound-free, free-free and bound-bound) opacity c fkhm = negative hydrogen ion opacity c fkm = molecular opacity c fkr = total radiative opacity c fkc = conductive opacity c fkap = total opacity (radiative and conductive) c -------- ad hoc correcting factors: A, B, C, introduced to beter reproduce c -------- results from more sophisticated program (with ionizations included) A=6.0 B=0.001 C=0.001 y=1-x-z zav=x+4*y+8*z fke=0.2*(1+x)/(1+2.7e11*rh/t/t)/(1+(t/4.5e8)**0.86) fkk=12.0/A*(1+x)*(0.001+z)*rh*(1.0e7/t)**3.5 fkhm=1.0e10 fkm=1.0e-5 if(t.gt.4.0e4)go to 1 fkhm=B*65.0*sqrt(z*rh)*(t/3000.0)**7.7 fkm=C*0.1*z 1 continue fkr=fkm+1.0/(1.0/fkhm+1.0/(fke+fkk)) fkc=1.0e10 if(rh.gt.0.00001) * fkc=zav*2.6e-7*(t/rh)**2*(1+(rh/2.0e6)**(2.0/3.0)) fkap=1.0/(1.0/fkr+1.0/fkc) return end ************* ******************** subroutine solve(deriv,delce,facdel) c ----- solves n=4 linear algebraic equations c INPUT: c deriv(n,n+1) = array of coefficients of "n" linear algebraic equations c OUTPUT: c delce(n) = n corrections, without delmax(i) they should satisfy n c equations: sum over j: deriv(i,j)*delce(j)+deriv(i,n+1)=0 c however, delce(i) are reduced so as to make them not larger c than delmax(i) c facdel = reduction factor for the corrections, if no reduction is c applied then facdel indicates the largest correction c - auxillary variables: c delmax(n) = maximum acceptable values of corrections c n = 4 = number of equations = number of unknowns (i.e. corrections) dimension deriv(4,5),delmax(4),delce(4) delmax(1)=0.03 delmax(2)=0.1 delmax(3)=0.01 delmax(4)=0.1 n=4 nm=n-1 np=n+1 do 1 k=1,nm kp=k+1 fac1=deriv(k,k) do 2 i=kp,n fac2=deriv(i,k) do 3 j=kp,np deriv(i,j)=deriv(i,j)*fac1-deriv(k,j)*fac2 3 continue 2 continue 1 continue c ------ the matrix of derivatives, deriv(i,j) is now triangular delce(n)=-deriv(n,np)/deriv(n,n) do 4 i=2,n i1=n-i+1 i2=i1+1 delce(i1)=-deriv(i1,np) do 5 j=i2,n delce(i1)=delce(i1)-deriv(i1,j)*delce(j) 5 continue delce(i1)=delce(i1)/deriv(i1,i1) 4 continue c ----- the unknowns delce(i) have been found c ----- check the size of "delce": dm=0.0 do 6 i=1,n d=abs(delce(i)/delmax(i)) if(dm.lt.d)dm=d 6 continue facdel=1.0 if(dm.gt.1.0)facdel=dm do 7 i=1,n delce(i)=delce(i)/facdel 7 continue facdel=dm return end ******************** ******************* subroutine prout(k,xx,ip,ic) c -- prints results of envelope or core integrations on the terminal and c -- on a unit number "ip"; if ic<0 then prints only the headings c INPUT: c k = integration step number, k<0 at the fitting point c xx(1) = temperature) (K) c xx(2) = density) (g/cm**3) c xx(3) = r (R Sun) c xx(4) = Lr (L Sun) c xx(5) = Mr (M Sun) c xx(6) = x = hydrogen content (by mass fraction) c xx(7) = z = heavy element content (by mass fraction) c xx(8) = zn = nitrogen content for CNO burning (by mass fraction) c ip = number of the output unit for "printing" c ic = control variable: ic>0 print results, ic<0 print table heading c ------ variables from common/const/ used by subroutine prout.for: c starm = total stellar mass / solar mass common/const/sunm,sunl,sunr,sunt,cgrrad,cdpm,sunmr3,sunml * ,pi,g,c,fmod dimension xx(10) if(ic.gt.0)go to 1 write(*,100) write(ip,100) 100 format(1x, * ' step log T log rho lg r/Rsun lg L/Lsun lg M/Msun') return 1 continue tlog=alog10(xx(1)) rhlog=alog10(xx(2)) rlog=alog10(xx(3)) fllog=alog10(xx(4)) fmlog=alog10(xx(5)) write(*,101)k,tlog,rhlog,rlog,fllog,fmlog write(ip,101)k,tlog,rhlog,rlog,fllog,fmlog 101 format(1x,i5,5f11.5) return end *****************