PROGRAM DWARF
C
C
C  ****************************************************************
C  THIS PROGRAM CALCULATES WHITE-DWARF STRUCTURES (FOR T = 0 K) FOR
C  ANY DEGREE OF RELATIVISTIC DEGENERACY.  AS INPUT IT REQUIRES THE
C  CENTRAL DENSITY IN G/CM^3  AND GIVES AS OUTPUT THE TOTAL MASS IN 
C  MSUN AND THE RADIUS OF THE OBJECT IN RSUN.
C  ****************************************************************
C
      IMPLICIT REAL*8 (A-H,O-Z)		
      DIMENSION YY(2)
C
C  DEFINITION OF CONSTANTS
C
      CMSUN=1.99E33  ! SOLAR MASS IN GRAM (NOT ACTUALLY USED)
      CRSUN=6.96E10  ! SOLAR RADIUS IN RSUN
      RHO0=3.789E6   ! MATCHING DENSITY IN G/CM^3
      ASCALE=1.557E8 ! LENGTH SCALE
C
 1    CONTINUE
      WRITE(*,*) ' CENTRAL DENSITY [G/CM^3]?'
      READ(*,*) RHOC
      THETAC=RHOC/RHO0  ! CONVERT INTO NON-DIMENSIONAL DENSITY
C
C  CENTRAL BOUNDARY CONDITIONS AND CONTROL PARAMETERS
C
      YY(1)=THETAC ! DENSITY VARIABLE
      YY(2)=0.     
      H=1.E-3 ! INITIAL STEP SIZE
      XEND=1.E-8  ! TO AVOID DIVISION BY ZERO AT CENTRE
      ACCC=1.E-10 ! ACCURACY PARAMETER (CONTROL PRECISION)
      ACCC0=1.E-11 ! HAS NO EFFECT IN THIS CASE
      N=2     ! NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS
      VMASS=0. ! CALCULATE NON-DIMENSIONAL MASS
C
C  ---------------------
C  MAIN INTEGRATION LOOP
C  ---------------------
C
 10   CONTINUE
      XX=XEND ! NEXT VALUE OF INDEPENDENT VARIABLE
      YY1O=YY(1)  ! STORE PREVIOUS VALUE OF DENSITY VARIABLE
      CALL MERS(XX,XEND,YY,N,ACCC,ACCC0,H) ! NEXT INTEGRATION STEP
      VMASS=VMASS+0.5*(XX*XX*YY1O+XEND*XEND*YY(1))*H ! MASS INTEGRATION 
      IF(YY(1)/THETAC.LT.0.001) GOTO 20 ! END INTEGRATION WHEN DENSITY IS LOW
      GOTO 10
C
 20   CONTINUE
      VR=XEND*ASCALE/CRSUN ! RADIUS IN SOLAR RADII
      VM=0.09*VMASS      ! MASS IN SOLAR MASSES
C
      WRITE(*,1000) ' RHOC = ',RHOC,', M/M_SUN = ',VM,', R/R_SUN = ',VR
 1000 FORMAT(A,1P,E10.4,0P,A,F7.4,A,F8.6)
      GOTO 1
C
      END