CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM DWARF C The structure of white dwarf stars CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INCLUDE 'DIM.P2' INCLUDE 'PAR.P2' C PRINT*,' enter value of central density (in units of gm/cm^3)' READ(*,*)RHO1 PRINT*,' enter number of electrons per nucleon' READ(*,*)YE C OPEN(12,FILE='LIXO.OUT', STATUS='UNKNOWN') CALL PARAM !get input parameters CALL ARCHON !calculate density, mass and radius of stars END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ARCHON C calculate the mass(r), density(r) for a given central density CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Local variables: REAL M !mass of star REAL R !radius of star REAL DR !radial step REAL EGRAV !gravitational energy of star REAL EKINT !kinetic energy of star REAL RHOCEN !central density INTEGER IR !number of radial steps C Functions: REAL GAMMA !dPressure/dDensity C Global variables: INCLUDE 'DIM.P2' REAL RHOSTR(0:MAXSTP) !density(r) of star INCLUDE 'PAR.P2' C C scaled central density RHOCEN=RHO1/RHO0 C DR=((3.0*0.001/RHOCEN)**(1.0/3.0))/3.0 !radial step R=DR/10. !begin at finite radius C C use Taylor expansion to obtain initial conditions M=(R**3.0)*(RHOCEN/3.0) !initial mass C initial density RHO=RHOCEN*(1-RHOCEN*R**2.0/6.0/GAMMA(RHOCEN**(1./3.))) RHOSTR(0)=RHO*RHO0 C C integrate equations and find energies, mass, and radius of star CALL INTGRT(M,RHO,R,DR,EGRAV,EKINT,IR,RHOSTR) C CALL TXTOUT(12,R,M,DR,EKINT,EGRAV,IR,RHOSTR) RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE INTGRT(M,RHO,R,DR,EGRAV,EKINT,IR,RHOSTR) C integrates mass and density begining with M and RHO at radius R; C calculates gravitational and kinetic energies CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Passed variables: REAL M,RHO !current mass and density (I/O) REAL R !current radius (I/O) REAL DR !radial step (input) REAL EGRAV !gravitational energy of star (output) REAL EKINT !kinetic energy of star (output) INTEGER IR !number of radial steps C Function: REAL EPSRHO !function in energy density C Global variables: INCLUDE 'DIM.P2' REAL RHOSTR(0:MAXSTP) !density(r) of star (output) INCLUDE 'PAR.P2' C EGRAV=0.0 !zero sums EKINT=0.0 IR=0 10 IF ((IR .LT. MAXSTP) .AND. (RHO .GT. RHOCUT)) THEN IR=IR+1 !update radial index CALL RUNGE(M,RHO,R,DR) !take a Runge-Kutta step IF (RHO .LT. RHOCUT) RHO=RHOCUT !avoid small densities RHOSTR(IR)=RHO*RHO0 !save values for output EGRAV=EGRAV+M*RHO*R !contribution to energy integ. EKINT=EKINT+R**2*EPSRHO(RHO**(1./3.)) GOTO 10 END IF EGRAV=-EGRAV*DR*E0 !unscaled energies EKINT=EKINT*DR*E0 RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE RUNGE(M,RHO,R,DR) C take a Runge-Kutta step to integrate mass and density CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Passed variables: REAL M,RHO !current mass and density (I/O) REAL R !current radius (input) REAL DR !radial step (input) C Local variables: REAL DMDR,DRHODR !radial derivatives of mass and dens REAL K1M,K1RHO !intermediate increments of mass and REAL K2M,K2RHO ! density REAL K3M,K3RHO REAL K4M,K4RHO C CALL DERS(M,RHO,R,DMDR,DRHODR) !calculate k1's K1M=DR*DMDR K1RHO=DR*DRHODR C R=R+DR/2.0 !calculate k2's M=M+.5*K1M RHO=RHO+.5*K1RHO CALL DERS(M,RHO,R,DMDR,DRHODR) K2M=DR*DMDR K2RHO=DR*DRHODR C M=M+.5*(K2M-K1M) !calculate k3's RHO=RHO+.5*(K2RHO-K1RHO) CALL DERS(M,RHO,R,DMDR,DRHODR) K3M=DR*DMDR K3RHO=DR*DRHODR C R=R+DR/2.0 !calculate k4's M=M+K3M-.5*K2M RHO=RHO+K3RHO-.5*K2RHO CALL DERS(M,RHO,R,DMDR,DRHODR) K4M=DR*DMDR K4RHO=DR*DRHODR C C values of mass and density at new radius M=M-K3M+(K1M+2.0*K2M+2.0*K3M+K4M)/6.0 RHO=RHO-K3RHO+(K1RHO+2.0*K2RHO+2.0*K3RHO+K4RHO)/6. RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE DERS(M,RHO,R,DMDR,DRHODR) C calculate the derivatives of mass and density with respect to radius CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Passed variables: REAL M,RHO !current mass and density (input) REAL R !current radius (input) REAL DMDR,DRHODR !radial derivatives of mass and dens (output) C Functions REAL GAMMA !dPressure/dDensity C IF (RHO .GT. 0.0) THEN !avoid division by zero DRHODR=-M*RHO/((R**2)*GAMMA(RHO**(1.0/3.0))) DMDR=(R**2)*RHO END IF RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION GAMMA(X) C derivative of pressure with respect to density CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Passed variables: REAL X !cubed root of density C GAMMA=X*X/3.0/SQRT(X*X+1) END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION EPSRHO(X) C function needed in the calculation of the energy CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Passed variables: REAL X !cubed root of density C Local variables: REAL PART1,PART2 !terms in the function C PART1=X*(1.0+2.0*X**2)*SQRT(1.0+X**2) PART2=LOG(X+SQRT(1.0+X**2)) EPSRHO=0.375*(PART1-PART2) END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE PARAM C calculate parameter dependent quantities CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global varibles: INCLUDE 'DIM.P2' INCLUDE 'PAR.P2' C R0=7.72E+08*YE !scaling radius M0=5.67E+33*(YE**2) !scaling mass RHO0=979000./YE !scaling density MEORMP=1./1836. E0=YE*MEORMP*9*(M0/1E+31) !scaling energy RHOCUT=1000./RHO0 !nearly zero density RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TXTOUT(MUNIT,R,M,DR,EKINT,EGRAV,IR,RHOSTR) C writes results to the requested unit CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Input variables: INTEGER MUNIT !unit to which we are writing REAL M !current mass and density REAL R !current radius REAL DR !radial step REAL EGRAV !gravitational energy of star REAL EKINT !kinetic energy of star INTEGER IR !number of radial steps C Global variables: INCLUDE 'DIM.P2' REAL RHOSTR(0:MAXSTP) !density(r) of star INCLUDE 'PAR.P2' C WRITE (MUNIT,19) WRITE (MUNIT,21) WRITE (MUNIT,27) YE WRITE (MUNIT,29) RHO1 WRITE (MUNIT,19) C WRITE (MUNIT,31) WRITE (MUNIT,32) WRITE (MUNIT,33) C C write unscaled values to output file WRITE (MUNIT,35) DR*R0,IR, + R*R0,M*M0,EKINT,EGRAV,EKINT+EGRAV C C write distance to center x density WRITE (MUNIT,19) WRITE (MUNIT,19) WRITE(MUNIT,36) WRITE(MUNIT,37) DO I=0,IR WRITE(MUNIT,38)I*DR*R0,RHOSTR(I) ENDDO C 19 FORMAT (' ') 21 FORMAT (' Output from project 2: ', + 'The structure of white dwarf stars ') 27 FORMAT (' Electron fraction YE= ',1PE9.3) 29 FORMAT (' Central density (gm/cm^3)=',1PE9.3) 31 FORMAT (3X,'Step',4X,'Nstep',3X,'Radius', + 5X,'Mass', 7X,'Kin E',6X,'Grav E',7X,'Tot E') 32 FORMAT (3X,' cm ',4X,' ',3X,' cm ', + 5X,' gm ', 7X,9X,'10^51 ergs ',4X,' ') 33 FORMAT (3X,'----',4X,'-----',3X,'------', + 5X,'----', 7X,'-----',6X,'------',7X,'-----') 35 FORMAT (1X,1PE8.2,1X,I4,3(2X,1PE9.3),2(2X,1PE10.3)) 36 FORMAT (10X,'Radius',12X,'Dens') 37 FORMAT (10X,'------',12X,'-----') 38 FORMAT (10X,1PE8.2,10X,1PE8.2) C RETURN END