program star2 ! ! to make use "f77 -g -n32 -L/usr/localn32/lib star2.f \ ! -lp329f -lodepack -llinpack -lblas -o star2 ! implicit none external fex, jdum integer neq parameter ( neq = 2 ) integer lrw, liw, jt, iout parameter (lrw = 22 + neq * 16) parameter (liw = 20 + neq) real*8 atol, rtol, rwork(lrw), r, rout, y(neq), tol real*8 C, lambda, mass, momen, dr, Pi, rho parameter (Pi = 3.141592654) common C, lambda real*8 r8arg integer iargc integer itol, iwork(liw), itask, istate, iopt, loopf ! ! Get arguments from the command line ! if (iargc() .eq. 5) then rho = r8arg(1,-1.0d0) C = r8arg(2,-1.0d0) lambda = r8arg(3, -1.0d0) dr = r8arg(4, -1.0d0) tol = r8arg(5, -1.0d0) else print *,"You must specify rho_0, C, lambda, dr and tol" stop endif ! !Check for invalid inputs ! if (rho .lt. 0) then print *,"rho_0 must be a positive number" stop endif if (dr .lt. 0) then print *,"dr must be a positive number" stop endif ! ! y(1) = m, y(2) = P ! y(1) = 0.000d0 y(2) = C*(rho**lambda) r = 0.0d0 rout = dr itol = 1 rtol = tol atol = tol itask = 1 istate = 1 iopt = 0 jt = 2 !integration loop 10 call lsoda(fex,neq,y,r,rout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jdum,jt) ! !Calculate rho from the pressure and print the results ! rho = (y(2)/C)**(1/lambda) write(6,30)r,rho,y(2),y(1) 30 format(4(e14.6, 2X)) ! ! Check for errors from lsoda ! if (istate .lt. 0) go to 80 rout = rout+dr ! !Stop if the pressure goes negative ! if (y(2) .gt. 0) then goto 10 else stop endif 60 format(2(e14.6, 2x)) 80 write(6,90)istate 90 format(///22h error halt.. istate =,i3) stop end ! !Subroutine for calculating the first derivatives ! subroutine fex (neq, r, y, ydot) implicit none integer neq double precision r, y(neq), ydot(neq), lambda, C, Pi, rho parameter (Pi = 3.141592654) common C, lambda ! !ydot(1) and ydot(2) are equal to dm/dr and dP/dr ! !Stop if pressure is less than zero ! if (y(2) .gt. 0) then ydot(1) = 4*Pi*(r**2)*((y(2)/C)**(1/lambda)) if (r .le. 0.0d0) then ydot(2) = 0.0d0 else ydot(2) = -(y(2)+((y(2)/C)**(1/lambda))) ydot(2) = ydot(2)*(y(1)+(4*Pi*(r**3)*y(2))) ydot(2) = ydot(2)/(r*(r-2.0d0*y(1))) endif ! !Write results to a file for checking ! open(1,'junk2') write (1,10)r,y(1),y(2),ydot(1),ydot(2) 10 format ('r= ',e14.6/ 1 'm= ',e14.6 ' P= ',e14.6/ 1 'mdot= ',e14.6 ' Pdot= ',e14.6/) return else rho = (y(2)/C)**(1/lambda) write(6,30)r,rho,y(2),y(1) 30 format(4(e14.6, 2X)) stop endif end ! ! Dummy jacobian routine: needed because there are no implicit variable decs ! subroutine jdum implicit none return end