{ Code © Charles Chandler http://qdl.scs-inc.us/?top=8583 } PROCEDURE Main; {This program uses finite element analysis to double-check the numbers in the Dalsgaard model. It is written to run as a plug-in for VectorWorks, but could be ported easily into ANSI Pascal.   The Sun is divided into 1000 pyramids, with "bases" equal to 1/1000 the surface area of the sphere, and which are then sub-divided into 100 truncated pyramids, which are then treated as discrete parcels in the FEA calcs. It doesn't matter that the square bases of the pyramids do not form a perfect polyhedron, because we just need to know the mass of each parcel, which we compress into the centroid. Then we can calculate the force of gravity, from each parcel, to each parcel. From this we get the pressure exerted by each parcel on each lower parcel. Knowing the pressure, and given the temperature, we can calculate the density of the 75/25 hydrogen/helium mix. The calculated densities are then compared to the initial values from the Dalsgaard model. The results are within 1/3 an order of magnitude. The nature of the discrepancy has not been determined. For one thing, I'm starting at 1.0 SR, while the Dalsgaard model has a bit of mass above the surface. The calcs also return an inordinate mass at the center. With more mass around the outside, and less in the center, the gravitational pressure would be less. Would it be 1/3 of an order of magnitude less? Only fixing the inaccuracies will tell for sure.   Note that to make the calcs run faster, it is not necessary to calculate the force of gravity from each parcel to each parcel. Rather, it is only necessary to calculate the forces acting on one pyramid, as all of them will be the same. So we take one pyramid and find the force of gravity from each of its parcels to each other parcel in each other pyramid. 1000 pyramids * 100 parcels = 100,000 calculations. But we don't need to repeat that process for the other 999 pyramids, totaling 99,900,000 calculations, because we'd get the same numbers every time.   Also note that this does not take the Coulomb barrier into account.   To see a graph of the results, * open 2ndParty/Images/Charles/Sun/Elements/Sun.vwx * click ReWriteGravities1 * copy Plug-ins/CLC/Geophysics/Data/1.Densities.txt to 2ndParty/Images/Charles/Sun/Dalsgaard/ * run the routine in Scratch.php that converts the results to logarithms * copy 1.Log(Densities).txt back to Plug-ins/CLC/Geophysics/Data/ * update the VSO in Sun.vwx } CONST {filenames} dataFileDensities0 = 'Plug-ins\CLC\Geophysics\Data\0.Log(Densities).txt'; {initial densities} dataFileTemperatures0 = 'Plug-ins\CLC\Geophysics\Data\0.Log(Temperatures).txt'; {initial temperatures} dataFileGravities1 = 'Plug-ins\CLC\Geophysics\Data\1.Gravities.txt'; {calculated gravities} dataFileDensities1 = 'Plug-ins\CLC\Geophysics\Data\1.Densities.txt'; {calculated densities} dataFileArchetype1 = 'Plug-ins\CLC\Geophysics\Data\1.Archetype.txt'; {calculated everything}   {ideal gas laws} uniGasConstant = 8.3144621; {J/K/mol} mole = 6.02214179 * (10^23); hydGramsPerAtom = 1.67262158 / (10^24); hydGramsPerMole = 1.00727643157738; helGramsPerAtom = 6.64648 / (10^24); helGramsPerMole = 4.00260449643992; mixGramsPerMole = (hydGramsPerMole * .75) + (helGramsPerMole * .25); hydGasConstant = 5621; {for plasma, extrapolated from other constants} helGasConstant = 2077; mixGasConstant = (hydGasConstant * .75) + (helGasConstant * .25);   {solar dimensions} solarRadius = 696; {radius of Sun, in Mm} solarVolume = 1412265429.1059; {for double-checking pyramids, in cubic Mm} solarMass = 1.9891 * (10^30); {kilograms} solarArea = 4 * 3.14159265 * solarRadius * solarRadius; {surface area of photosphere, in square Mm}   {programmatic constants} dataLineCnt = 2482; numPyramids = 1000; {number of pyramidal subdivisions of Sun} pyramidBaseArea = solarArea / numPyramids; {area of each pyramid base, in square Mm} pyramidBaseWidth = SqRt(pyramidBaseArea); {width of pyramid, as square root of area} pyramidVolume = (1/3) * pyramidBaseArea * solarRadius; {volume of each pyramid, in cubic Mm} numPyramidSlices = 100; {number of truncated pyramids in each stack} aveDensityError100 = 1362.5379796089; {at 100 pyramid slices, the average density is 96.77% of 1408} aveDensityError1000 = 1401.68838482822; {at 1000 pyramid slices, the average density is 99.55% of 1408} pyramidSliceVolume = pyramidVolume / numPyramidSlices; {volume of each truncated pyramid, in cubic Mm} pyramidBaseOverHght = (pyramidBaseWidth / 2) / solarRadius; {opposite / adjacent, faster than trig} pyramidAngle = ArcTan(1 / pyramidBaseOverHght); {slope of pyramid from horizontal, in radians} TYPE parcel = STRUCTURE outerRad :REAL; {outer radius} centroid :REAL; {radius at the centroid} baseArea :REAL; {area of truncated pyramid base} volume :REAL; {volume of truncated pyramid} density :REAL; {initialized from Dalsgaard model, then recalculated} mass :REAL; {volume * density, used in gravity calcs} gravity :VECTOR; {sum of gravity vectors} pressure :REAL; {gravity per lower area} temp :REAL; {temperature} atomNum :REAL; {atomic number of element} END; dataLine = STRUCTURE radius :REAL; {r/R} speed :REAL; {cm/sec} density :REAL; {g/cm^3} pressure :REAL; {dyn/cm^2} gamma :REAL; {adiabats} temp :REAL; {K} END; VAR objName :STRING; objHand, recHand, wallHand :HANDLE; spherePts :ARRAY [1..numPyramids] OF VECTOR; onePyramid :ARRAY [1..numPyramidSlices] OF parcel; dataLines :ARRAY [1..dataLineCnt] OF dataLine; dataLinesCursor :INTEGER; gramsPerMole :ARRAY [1..118] OF REAL;   h :HANDLE; str :STRING; inc, off, x, y, z, r, phi, tmpReal :REAL; cnt, cnt1, cnt2, cnt3 :INTEGER; runningVolume, runningGravity :REAL; baseArea, baseWidth, peakArea, height, hypeScaleX, hypeScaleY :REAL; one, archetype1, archetype2 :parcel; pt1, pt2, tmpVect :VECTOR; tmpNorm, tmpGrav, tmpMoles :REAL; tick1, tick2 :LONGINT;   FUNCTION ForceOfGravity(m1, m2, dist :REAL) :REAL; BEGIN {Given the mass (kg) and distance (m), return the gravitational force (N). http://easycalculation.com/physics/classical-physics/newtons-law.php} ForceOfGravity := ((6.6726 / (10^11)) * m1 * m2) / (dist^2); END;   PROCEDURE ForceOfGravityTest; VAR tmpReal :REAL; BEGIN {These are numbers from somewhere else, used to double-check the gravity calcs.} tmpReal := ForceOfGravity(222222, 222222, 1); IF (Abs(tmpReal - 3.2951045208922) > .00001) THEN Message('error: gravity calculator discrepancy: ', tmpReal); tmpReal := ForceOfGravity(2222222, 2222222, 2); IF (Abs(tmpReal - 82.377761302223) > .00001) THEN Message('error: gravity calculator discrepancy: ', tmpReal); END;   FUNCTION HeightOfPyramid(ang, vol: REAL) :REAL; BEGIN {http://www.cleavebooks.co.uk/scol/calpyr.htm} {ang is measured from the horizontal.} HeightOfPyramid := ((6 * vol / Tan(ang)) ^ (1/3)) * Tan(ang) / 2; END;   FUNCTION VolumeOfPyramid(a, b, h: REAL) :REAL; BEGIN {http://mathworld.wolfram.com/TruncatedSquarePyramid.html} {a = bottom, b = top, h = height} {For a non-truncated pyramid, b = 0, and it's just 1/3 * a * a * h.} VolumeOfPyramid := (1/3) * ((a * a) + (a * b) + (b * b)) * h; END;   FUNCTION CentroidOfPyramid(a1, a2, h :REAL) :REAL; BEGIN {http://mathworld.wolfram.com/PyramidalFrustum.html} {a1 = bottom area, a2 = top area, h = height} CentroidOfPyramid := (h * (a1 + (2 * SqRt(a1 * a2)) + (3 * a2))) / (4 * (a1 + SqRt(a1 * a2) + a2)); END;   PROCEDURE PlaceText(text :STRING; x, y, rot, size :REAL; just, vAlign :STRING); CONST isMirrored = false; VAR justInt, vAlignInt :INTEGER; BEGIN CreateText(text); SetTextOrientation(LNewObj, x, y, rot, isMirrored); IF (just = 'left' ) THEN justInt := 1 ELSE IF (just = 'center') THEN justInt := 2 ELSE IF (just = 'right' ) THEN justInt := 3; SetTextJust(LNewObj, justInt); IF (vAlign = 'top' ) THEN vAlignInt := 1 ELSE IF (vAlign = 'textTop') THEN vAlignInt := 2 ELSE IF (vAlign = 'center' ) THEN vAlignInt := 3 ELSE IF (vAlign = 'textBot') THEN vAlignInt := 4 ELSE IF (vAlign = 'bot' ) THEN vAlignInt := 5; SetTextVerticalAlign(LNewObj, vAlignInt); SetTextSize(LNewObj, 0, Len(text), size); END;   FUNCTION LPad(theString :STRING; theLength :INTEGER; theCharacter :STRING) :STRING; BEGIN WHILE Len(theString) < theLength DO theString := Concat(theCharacter, theString); LPad := theString; END;   BEGIN IF GetCustomObjectInfo(objName, objHand, recHand, wallHand) THEN BEGIN {Double-check total volume.} ClrMessage; tmpReal := (pyramidVolume * numPyramids) - solarVolume; IF (Abs(tmpReal) > 1.6138) THEN Message('error: overall volume discrepancy: ', tmpReal);   {Get the 3D points on the surface of the sphere, using the Golden Spiral method. Note that this misses the last one, so we tell it to create (numPyramids - 1), then we create the last one manually, bringing the total up to numPyramids.} inc := PI * (3 - sqrt(5)); off := 2 / (numPyramids - 1); FOR cnt := 1 TO (numPyramids - 1) DO BEGIN y := (cnt * off) - 1 + (off / 2); r := SqRt(Abs(1 - (y * y))); phi := cnt * inc; spherePts[cnt].x := cos(phi) * r; spherePts[cnt].y := y; spherePts[cnt].z := sin(phi) * r; END; {Throw in the one that it misses.} spherePts[numPyramids].x := .05; spherePts[numPyramids].y := -1; spherePts[numPyramids].z := 0; {Make sure the basePts are exactly unit vectors (length = 1).} FOR cnt := 1 TO numPyramids DO spherePts[cnt] := UnitVec(spherePts[cnt]);   {Load the log(densities) from the text file.} Open(dataFileDensities0); FOR cnt := dataLineCnt DOWNTO 1 DO BEGIN ReadLn(dataLines[cnt].radius, dataLines[cnt].density); dataLines[cnt].density := (10 ^ dataLines[cnt].density) * 1000; {convert g/cm^3 to kg/m^3} END; Close(dataFileDensities0);   {Load the log(temperatures) from the text file.} Open(dataFileTemperatures0); FOR cnt := dataLineCnt DOWNTO 1 DO BEGIN ReadLn(dataLines[cnt].radius, dataLines[cnt].temp); dataLines[cnt].temp := (10 ^ dataLines[cnt].temp); END; Close(dataFileTemperatures0);   {Assign the info for the archetype pyramid.} dataLinesCursor := 1; FOR cnt := 1 TO numPyramidSlices DO BEGIN IF (cnt = 1) THEN BEGIN {Create a tiny parcel at the very bottom, just so we get accurate numbers there, not just extrapolations of the innermost parcel, which is very skinny & tall.} onePyramid[cnt].outerRad := 1; onePyramid[cnt].volume := VolumeOfPyramid(pyramidBaseOverHght * 2, 0, 1); END ELSE IF (cnt = 2) THEN BEGIN {The second one is a standard slice, but with the innermost parcel subtracted from it.} runningVolume := pyramidSliceVolume; onePyramid[cnt].outerRad := HeightOfPyramid(pyramidAngle, runningVolume); onePyramid[cnt].volume := pyramidSliceVolume - onePyramid[1].volume; END ELSE BEGIN runningVolume := runningVolume + pyramidSliceVolume; onePyramid[cnt].outerRad := HeightOfPyramid(pyramidAngle, runningVolume); onePyramid[cnt].volume := pyramidSliceVolume; END; baseWidth := onePyramid[cnt].outerRad * pyramidBaseOverHght * 2; onePyramid[cnt].baseArea := baseWidth * baseWidth;   {Find the centroid.} baseArea := onePyramid[cnt].baseArea; IF (cnt = 1) THEN BEGIN peakArea := 0; height := onePyramid[cnt].outerRad; END ELSE BEGIN peakArea := onePyramid[cnt-1].baseArea; height := onePyramid[cnt].outerRad - onePyramid[cnt-1].outerRad; END; onePyramid[cnt].centroid := onePyramid[cnt].outerRad - CentroidOfPyramid(baseArea, peakArea, height);   {Find the density at the centroid and calculate the mass.} tmpReal := onePyramid[cnt].centroid / solarRadius; {find the r/R for this slice} WHILE (tmpReal > dataLines[dataLinesCursor].radius) DO dataLinesCursor := dataLinesCursor + 1; onePyramid[cnt].mass := dataLines[dataLinesCursor - 1].density * onePyramid[cnt].volume * (10^18); onePyramid[cnt].temp := dataLines[dataLinesCursor - 1].temp;   {Just initialize these.} onePyramid[cnt].gravity.x := 0; onePyramid[cnt].gravity.y := 0; onePyramid[cnt].gravity.z := 0; onePyramid[cnt].pressure := 0; onePyramid[cnt].density := 0; END;   {Double-check the mass. It actually comes out lower than the nominal 1408 kg/m^3, perhaps 'cuz of the way I'm applying the centroid density to the whole parcel. But within 97% is good enough.} tmpReal := 0; FOR cnt := 1 TO numPyramidSlices DO tmpReal := tmpReal + onePyramid[cnt].mass; tmpReal := (tmpReal * 1000) / (solarVolume * (10^18)); {get the kg/m^3} IF ((numPyramidSlices = 100 ) AND (Abs(tmpReal - aveDensityError100 ) > 1)) THEN Message('error: overall density: ', tmpReal); IF ((numPyramidSlices = 1000) AND (Abs(tmpReal - aveDensityError1000) > 1)) THEN Message('error: overall density: ', tmpReal);   IF pReWriteGravities1 THEN BEGIN {Find the force of gravity, from each parcel, to each parcel, adding the vector to what's already there. For the archetypical parcels, we use spherePts[500], because it is in the middle of the generated points, where the points are evenly spaced, as opposed to the beginning or end of the Golden Spiral, where things get a little irregular, and nearby irregularities would throw off the gravity calcs. At 1000 slices, after 2 hours, it had only done 416 slices. At 100 slices, seconds elapsed: 161} tick1 := GetTickCount; ReWrite(dataFileGravities1); FOR cnt1 := 1 TO numPyramidSlices DO BEGIN archetype1 := onePyramid[cnt1]; IF (archetype1.gravity.x = 0) THEN BEGIN pt1 := spherePts[500] * archetype1.centroid * (10^6); FOR cnt2 := 1 TO numPyramidSlices DO BEGIN archetype2 := onePyramid[cnt2]; FOR cnt3 := 1 TO numPyramids DO BEGIN {Calculate the gravity from each slice to each slice, in each of the pyramids, but skip it if it's actually the same parcel, saving 100 calcs.} IF ((cnt1 <> cnt2) OR (cnt3 <> 500)) THEN BEGIN pt2 := spherePts[cnt3] * archetype2.centroid * (10^6); tmpVect := pt2 - pt1; {vector between centroids} tmpNorm := Norm(tmpVect); {magnitude of that vector} tmpGrav := ForceOfGravity(archetype1.mass, archetype2.mass, tmpNorm); tmpVect := UnitVec(tmpVect) * tmpGrav; {gravity as a vector} onePyramid[cnt1].gravity := onePyramid[cnt1].gravity + tmpVect; END; END; END; END; WriteLn(Concat( onePyramid[cnt1].gravity.x, Chr(9), onePyramid[cnt1].gravity.y, Chr(9), onePyramid[cnt1].gravity.z)); END; Close(dataFileGravities1); tick2 := GetTickCount; WriteLn('seconds elapsed: ', (tick2 - tick1) / 60); SetRField(objHand, objName, 'ReWriteGravities1', 'False'); END ELSE BEGIN {To speed up debugging, we can skip the calculation of the gravity vectors, which takes several minutes, by just re-loading values previously generated and stored in a text file, assuming the gravity calcs don't need debugging.} dataLinesCursor := 0; Open(dataFileGravities1); WHILE NOT EOF(dataFileGravities1) DO BEGIN ReadLn(x, y, z); IF (Abs(x) > 0) THEN BEGIN dataLinesCursor := dataLinesCursor + 1; onePyramid[dataLinesCursor].gravity.x := x; onePyramid[dataLinesCursor].gravity.y := y; onePyramid[dataLinesCursor].gravity.z := z; END; END; Close(dataFileGravities1); END;   {Check that the gravity vector is pointing at the center of gravity. The innermost parcel is off (0.14), but its gravity doesn't matter.} FOR cnt := 2 TO numPyramidSlices DO BEGIN tmpReal := Norm(UnitVec(onePyramid[cnt].gravity) + spherePts[500]); IF (tmpReal > 0.002) THEN Message('error: gravity vector: ', cnt); END;   {Calculate the pressure, given the force of gravity acting on the parcels above it. Note that gravity is stored as a vector, which is the sum of all of the centroid-to-centroid forces. This always ends up pointing at the center of gravity, due to the spherical geometry. So in the end, we just need to know the magnitude of the vector, which is the gravitational force that parcel exerts on the parcel below it. The magnitude is in Newtons, while pressure is in N/m^2, so we divide by the base area, which is stored in square Mm.} IF (true) THEN BEGIN {This method cumulatively adds the pressures. The top parcel doesn't have anything above it, so it is held down just by the force of gravity. In actuality, only the topmost atoms have just the force of gravity holding them down, while at the bottom of the parcel, this pressure has accumulated. This will be true no matter how small the parcel (above the atomic level). So we just set the first parcel equal to 1.} onePyramid[numPyramidSlices].pressure := 1; FOR cnt := numPyramidSlices - 1 DOWNTO 1 DO BEGIN onePyramid[cnt].pressure := (Norm(onePyramid[cnt + 1].gravity) / (onePyramid[cnt].baseArea * (10^12))) + onePyramid[cnt + 1].pressure; END; END ELSE BEGIN {This way cumulatively adds the gravities. This produces slightly larger values.} runningGravity := Norm(onePyramid[numPyramidSlices].gravity); FOR cnt := numPyramidSlices - 1 DOWNTO 1 DO BEGIN onePyramid[cnt].pressure := (runningGravity / (onePyramid[cnt].baseArea * (10^12))); runningGravity := runningGravity + Norm(onePyramid[cnt].gravity); END; END;   {Knowing the pressures and temperatures, calculate the densities.} ReWrite(dataFileDensities1); FOR cnt := 1 TO numPyramidSlices DO BEGIN {ideal gas law: PV = nRT moles = (pressure * volume) / (gasConstant * temperature) massKg = (moles * mixGramsPerMole) / 1000 density = massKg / volume} tmpMoles := (onePyramid[cnt].pressure * onePyramid[cnt].volume) / (uniGasConstant * onePyramid[cnt].temp); onePyramid[cnt].mass := (tmpMoles * mixGramsPerMole) / 1000; {kg} onePyramid[cnt].density := onePyramid[cnt].mass / onePyramid[cnt].volume;   {Do it the easy way, just to double-check. https://www.brisbanehotairballooning.com.au/faqs/education/116-calculate-air-density.html} tmpReal := onePyramid[cnt].pressure / (onePyramid[cnt].temp * mixGasConstant); IF (Abs(1 - (tmpReal / onePyramid[cnt].density)) > 0.0001) THEN Message('error: ideal gas law double-check');   WriteLn(Concat( onePyramid[cnt].outerRad / solarRadius, Chr(9), onePyramid[cnt].density / 1000)); {convert kg/m^3 to g/cm^3} END; Close(dataFileDensities1);   {Write onePyramid to a text file.} ReWrite(dataFileArchetype1); WriteLn(Concat( LPad('radius (r/R)', 12, ' '), LPad('base (Mm^2)', 16, ' '), LPad('volume (Mm^3)', 16, ' '), LPad('density (kg/m^3)', 32, ' '), LPad('mass (kg)', 32, ' '), LPad('gravity (N)', 32, ' '), LPad('pressure (N/m^2)', 32, ' ') )); FOR cnt := numPyramidSlices DOWNTO 1 DO BEGIN one := onePyramid[cnt]; WriteLn(Concat( LPad(Num2Str(8, one.outerRad / solarRadius), 12, ' '), LPad(Num2Str(8, one.baseArea), 16, ' '), LPad(Num2Str(3, one.volume), 16, ' '), LPad(Num2Str(8, one.density), 32, ' '), LPad(Num2Str(0, one.mass), 32, ' '), LPad(Num2Str(0, Norm(one.gravity)), 32, ' '), LPad(Num2Str(8, one.pressure), 32, ' ') )); END; Close(dataFileArchetype1);   {Generate display graphics.} IF (pDisplay = 'Cones') THEN BEGIN {Display the pyramids as cones, to visually check the spacing.} Locus3D(spherePts[1].x, spherePts[1].y, spherePts[1].z); FOR cnt := 1 TO numPyramids DO BEGIN h := CreateCone(spherePts[cnt].x, spherePts[cnt].y, spherePts[cnt].z, 0, 0, 0, .05); END; END ELSE IF (pDisplay = 'Slices') THEN BEGIN FOR cnt := 1 TO numPyramidSlices DO BEGIN {Locus(0, onePyramid[cnt].centroid);} str := Concat(' ', Num2Str(4, onePyramid[cnt].density), ' kg/m^3, ', onePyramid[cnt].mass, ' kg'); {PlaceText(str, 0, onePyramid[cnt].centroid, 0, 3, 'left', 'center');} y := onePyramid[cnt].outerRad; x := y * pyramidBaseOverHght; IF (cnt MOD 5 = 0) THEN MoveTo(-x, y); LineTo(x, y); END; LineTo(0, 0); LineTo(-x, y); END ELSE IF (pDisplay = 'Spiral') THEN BEGIN {h := CreateNurbsCurve(spherePts[1].x, spherePts[1].y, spherePts[1].z, true, 1); FOR cnt := 2 TO numPyramids DO BEGIN AddVertex3D(h, spherePts[cnt].x, spherePts[cnt].y, spherePts[cnt].z); END;}   FOR cnt := 1 TO numPyramids DO BEGIN Locus3D(spherePts[cnt].x, spherePts[cnt].y, spherePts[cnt].z); END;   {BeginPoly3D; FOR cnt := 1 TO numPyramids DO BEGIN IF (cnt MOD 2 = 0) THEN Add3DPt(spherePts[cnt].x, spherePts[cnt].y, spherePts[cnt].z); END; EndPoly3D; SetFPat(LNewObj, 0);} END; END; END; RUN(Main);