program Planet (input*, output); const n = 2; pi = 3.1415926; beeping = true; type index = 1..n; yesno = (yes, no); vector = record c: array [index] of real end; pvec = ^ vector; plotrec = record maxx, minx, maxy, miny: real end; pplot = ^ plotrec; var trac: yesno; plota: text; plotb: text; ppa, ppb: pplot; procedure Beep (c: char); begin if beeping then write (output, c, ' '); end; procedure NewCurve (pp: pplot; var f: text); begin writeln(f); writeln (f, 'New curve'); writeln (f, 'Curve Type Solid'); writeln (f, 'Curve Interpolation Linear'); writeln (f, 'Add ? points'); end; procedure PlotAxis (var f: text; ax: char; mn, mx: real); var un,int: real; nd: integer; begin int := (mx-mn)/15.0; if int <= 0.0 then int := 1.0 else begin nd := round(ln(int)/ln(10.0) - 0.5); if nd = 0 then un := 1.0 else un := exp (nd * ln(10.0)); writeln(output, 'un===',un); break(output); int := trunc (int/un); if int > 5.0 then begin un := 10.0*un; int := un end else int := int*un; end; writeln(f, ax, ' minimum ', mn); writeln(f, ax, ' maximum ', mx+int); writeln(f, ax, ' interval ', int); write(f, ax, ' numberstyle ', max(0, -nd)); if nd < 0 then writeln(f, ',F') else writeln(f); end; procedure PlotPoint (pp: pplot; var f: text; x,y: real); begin writeln (f, x, ', ', y); pp^.maxx := max (pp^.maxx, x); pp^.minx := min (pp^.minx, x); pp^.maxy := max (pp^.maxy, y); pp^.miny := min (pp^.miny, y) end; procedure OpenPlot (pp: pplot; var f: text); begin rewrite (f); writeln(f, 'New Graph'); writeln(f, 'Size 5 by 8'); writeln(f, 'Orientation Landscape'); writeln(f, 'Graph font HELVETICA12BI'); writeln(f); pp^.maxx := -1000.0; pp^.minx := 1000.0; pp^.maxy := -1000.0; pp^.miny := 1000.0; end; procedure ClosePlot (pp: pplot; var f: text); begin writeln(f); PlotAxis (f, 'x', pp^.minx, pp^.maxx); PlotAxis (f, 'y', pp^.miny, pp^.maxy); writeln (f, 'plot dover,plot.press,print'); end; procedure Integrate ( function F(real; pvec; pvec): boolean; (* computes derivative *) xin: real; (* initial x *) y: pvec; (* initial/final function value *) h: real; (* integration step *) var xfin: real); (* final x *) (* (* Integrates the differential equation y = F (x, y) by the Runge-Cooka (* method. (* (* The integration starts at x = xin with the given initial value for y. (* It stops when x = xfin, or when F returns true;in the latter case, (* xfin is set to the last x for which y was computed. (* (* The function F should take x and y as the first two parameters, (* and return dy/dx in the third. Instability, stiffness, inaccuracy (* and similar problems are solved by the Ostrich method. (* *) label 9; var yhf, ydot, ydothf: pvec; i: index; nit: integer; hfh, x, yhhi, yi: real; stop: boolean; procedure TracePoint (nx: integer; x: real; y, yd: pvec); var j: index; begin writeln (output); write(output, '[INTEG] X', nx:3, ': ', x:16:8, ' Y: '); for j := 1 to n do begin write (output, y^.c[j] :16:8) end; writeln (output); write(output, ' ', ' D: '); for j := 1 to n do begin write (output, yd^.c[j] :16:8) end end; (* TracePoint *) begin (* Integrate *) new(yhf); new(ydot); new(ydothf); x := xin; hfh := h/2.0; nit := 0; while x+h <= xfin do begin stop := F (x, y, ydot); nit := nit + 1; if stop then goto 9; if trac=yes then TracePoint (nit, x, y, ydot); for i := 1 to n do yhf^.c[i] := y^.c[i] + hfh * ydot^.c[i]; stop := F (x + hfh, yhf, ydothf); nit := nit + 1; if stop then begin x := x + hfh; for i := 1 to n do y^.c[i] := yhf^.c[i]; goto 9 end; if trac=yes then TracePoint (nit, x, yhf, ydothf); for i := 1 to n do begin yhhi := yhf^.c[i] + ydothf^.c[i] * hfh; yi := y^.c[i] + ydot^.c[i] * h; y^.c[i] := 2.0 * yhhi - yi end; x := x + h; end; if x <= xfin then begin stop := F (x, y, ydot); nit := nit + 1; if stop then goto 9; if trac=yes then TracePoint (nit, x, y, ydot); for i := 1 to n do y^.c[i] := y^.c[i] + (xfin - x) * ydot^.c[i]; x := xfin end; 9: xfin := x end; (* Integrate *) procedure InteTest; var xin, xfin, h: real; i: integer; y: pvec; function Sinfun (x: real; y, ydot: pvec): boolean; (* a sinusoidal function *) begin Sinfun := false; ydot^.c[1] := y^.c[2]; ydot^.c[2] := -y^.c[1] end; (* NFunc *) procedure SinfunInit (var xin, xfin, h: real; y: pvec); begin xin := 0.0; xfin := pi; h := pi/10.0; y^.c[1] := 0.0; y^.c[2] := 1.0; end; begin (* InteTest *) new(y); writeln(output, ' FULL TRACE (YES OR NO)?'); break(output); read(input, trac); SinFunInit (xin, xfin, h, y); Integrate (SinFun, xin, y, h, xfin); writeln (output); write (output, '[INTEG] XFIN :', xfin:16:8, ' Y: '); for i := 1 to n do write (output, y^.c[i]:16:8); writeln (output) end; (* InteTest *) procedure LaneEmden; (* Integrates the Lane-Emden equations for a polytropic gas sphere *) const maxnkap = 10; var kaps: array [1..maxnkap] of real; kap, kap1: real; (* polytropic coefficient *) xin, xfin, h: real; y: pvec; i, j, nkap: integer; function LanEm (x: real; y, ydot: pvec): boolean; (* the Lane-Emden equations *) var p, q, w, s, u, pdot, qdot, ex: real; begin LanEm := false; s := x; p := y^.c[1]; q := y^.c[2]; if trac=yes then writeln (output, 'S=', s, 'P=', p, 'Q=', q); if x = 0.0 then begin pdot := 2.0*kap1; qdot := 1.0/(pdot*sqr(pdot)); u := 0.0; w := 0.0; end else begin if kap <= 0.0 then ex := 1.0 else ex := exp(-kap*s); pdot := 2.0*kap1*sqrt(p/s*sqr(p/s)*q/s)*ex; qdot := (q/s)*(4.0 - 6.0*kap1*sqr(p/s)*q/s*exp(-s)*ex); u := sqrt(p); w := sqr(s)/sqrt(q); end; if trac=yes then writeln (output, 'DP=', pdot, 'DQ=', qdot); ydot^.c[1] := pdot; ydot^.c[2] := qdot; PlotPoint (ppa, plota, u, s/ln(10.0)); PlotPoint (ppb, plotb, u, w); end; (* NFunc *) procedure LanEmInit (var xin, xfin, h: real; y: pvec); begin write (output, ' KAPPA = '); break(output); read (input, kap); kap1 := kap+1.0; write (output, ' MAX S (BASE 10) = '); break(output); read(input, xfin); xfin := xfin * ln(10.0); NewCurve (ppa, plota); NewCurve (ppb, plotb); xin := 0.0; h := ln(10.0)/20.0; y^.c[1] := 0.0; y^.c[2] := 0.0; end; procedure LanEmFinis; var j: integer; begin end; begin (* LaneEmden *) new(y); OpenPlot(ppa, plota); OpenPlot(ppb, plotb); write(output, ' N CURVES = '); break (output); read (input, nkap); for j := 1 to nkap do begin LanEmInit (xin, xfin, h, y); kaps[j] := kap; Integrate (LanEm, xin, y, h, xfin); writeln (output); write (output, '[INTEG] XFIN :', xfin:16:8, ' Y: '); for i := 1 to n do write (output, y^.c[i]:16:8); writeln (output) end; write(plota, 'Graph Title Co-log density for K ='); for j := 1 to nkap do write(plota, kaps[j]:7:3, ' '); writeln(plota); write(plotb, 'Graph Title Mass fraction for K ='); for j := 1 to nkap do write(plotb, kaps[j]:7:3, ' '); writeln(plotb); writeln(plota, 'X Label Normalized radius'); writeln(plotb, 'X Label Normalized radius'); ClosePlot (ppa, plota); ClosePlot (ppb, plotb); end; (* InteTest *) begin (* MAIN *) new(ppa); new(ppb); write (output, 'Trace (yes or no)?'); break(output); read (input, trac); LaneEmden; end.