DEFDBL a-z DECLARE FUNCTION f (x) DECLARE FUNCTION lg10 (x) DECLARE FUNCTION dfdx (x) DECLARE SUB density (p, xold, rho, xnew) DECLARE FUNCTION press (m, r, d) DECLARE FUNCTION mass (d, r) COMMON SHARED g, m0, r0, a, pi ' some physical quantities used in the program g = 6.673E-08 m0 = 2E+33 r0 = 6.96E+10 a = 6.01E+22 pi = 3.1415926536# CLS PRINT "Astrophysics with a PC : WHITE DWARF" PRINT "------------------------------------" PRINT " " PRINT "----- Minimal solution program -----" PRINT " " PRINT "Input of initial conditions and model parameters : " INPUT "Log10 of the central density : ", rhc INPUT "stepsize in km : ", drkm ' central values of the physical variables rhoc = 10 ^ rhc dr = 100000! * drkm PRINT " " m = 0! r = 0! xc = (rhoc / 1964000!) ^ (1 / 3) pc = a * f(xc) i% = 0 ' display heading of the table with the results PRINT " i r Mr log(P) log(rho) x" ' display central values (layer zero) PRINT USING "### ###.##### ###.##### ###.##### " &_ "###.##### ##.#######";_ i%; r / r0; m / m0; lg10(pc); lg10(rhoc); xc ' compute first step from layer zero to layer one ' and show them on screen p = pc - 2 / 3 * g * pi * (rhoc * dr) ^ 2 m = 4 / 3 * pi * rhoc * dr ^ 3 CALL density(p, xc, d, x) r = dr i% = 1 PRINT USING "### ###.##### ###.##### ###.##### " &_ "###.##### ##.#######";_ i%; r / r0; m / m0; lg10(p); lg10(d); x ' prepare for the other layers(from layer two to the surface) surface% = 0 i% = 2 ' here starts the main cycle for the other layers DO r1 = r + .5 * dr p1 = p + .5 * dr * press(m, r, d) ' check if pressure (half step i+1/2) is still positive ' otherwise surface is reached IF p1 > 0 THEN m1 = m + .5 * dr * mass(d, r) CALL density(p1, x, d1, x1) r = r + dr p = p + dr * press(m1, r1, d1) ' check if pressure (full step i+1) is still positive ' otherwise surface is reached IF p > 0 THEN m = m + dr * mass(d1, r1) CALL density(p, x1, d, x) ' show results of newly computed layer on screen PRINT USING "### ###.##### ###.##### " &_ "###.##### ###.##### ##.#######"; _ i%; r / r0; m / m0; lg10(p); lg10(d); x i% = i% + 1 ' pause after every 12 layers IF i% MOD 12 = 0 THEN PRINT "Press any key to continue" DO LOOP WHILE INKEY$ = "" END IF ELSE 'this ELSE is reached if pressure (i+1) 'was negative surface% = 1 END IF ELSE 'this ELSE is reached if pressure (i+1/2) 'was negative surface% = 1 END IF LOOP UNTIL surface% = 1 ' show general results on screen PRINT " " PRINT USING "Total mass : ##.##"; m / m0 PRINT "Program terminated. Press any key to stop" DO LOOP WHILE INKEY$ = "" END SUB density (p, xold, rho, xnew) ' ' computes the density rho and the corresponding ' value xnew from the ' input pressure P and starting guess xold k = p / a stopcrit% = 0 DO xnew = xold - (f(xold) - k) / dfdx(xold) IF ABS(xold - xnew) < .000001 THEN stopcrit% = 1 ELSE xold = xnew END IF LOOP UNTIL stopcrit% = 1 rho = 1964000! * xnew ^ 3 END SUB FUNCTION dfdx (x) ' ' evaluates the derivative of f(x) ' dfdx = 8 * x ^ 4 / SQR(1 + x ^ 2) END FUNCTION FUNCTION f (x) ' ' evaluates the function f(x) ' f = x * (2 * x ^ 2 - 3) * SQR(1 + x ^ 2) + _ 3 * LOG(x + SQR(1 + x ^ 2)) END FUNCTION FUNCTION lg10 (x) ' ' computes the base 10 logarithm of x using the natural ' logarithm LOG(x) ' lg10 = LOG(x) / LOG(10!) END FUNCTION FUNCTION mass (d, r) ' ' computes the right hand side of the differential equation ' of mass continuity ' mass = 4 * pi * d * r * r END FUNCTION FUNCTION press (m, r, d) ' ' computes the right hand side of the equation of hydro- ' static equilibrium ' press = -g * m * d / r ^ 2 END FUNCTION