{ New Code © Charles Chandler http://qdl.scs-inc.us/?top=9719 } PROCEDURE Main; CONST dataFile = 'Plug-ins\CLC\Geophysics\Data\Solar Elements.txt'; debugOutFile = 'debug.txt'; debug = false; lastIndex = 82; xMax = lastIndex * 2; bezPrec = lastIndex; TYPE row = STRUCTURE sym :STRING; {element abbr} nam :STRING; {element name} den :REAL; {density} apv :REAL; {atoms per volume} set1 :REAL; {Photosphere, from Anders & Grevesse (1989)} set2 :REAL; {Solar Wind, from Manuel, Kamat, & Mozina (2006)} set3 :REAL; {B2FH Values, from Manuel, Kamat, & Mozina (2006)} use :REAL; {which set to use} cor :REAL; {correction} vol :REAL; {volume} rad :REAL; {radius} END; VAR objName, noneClass, str :STRING; objHand, recHand, wallHand, pathHand, correctPoly :HANDLE; rows :ARRAY [0..lastIndex] OF row; bezier :ARRAY [1..bezPrec] OF POINT; p0, p1, p2, p3 :POINT; cnt, cnt1, cnt2, vertexType, zoomFactor :INTEGER; p1x, p1y, p2x, p2y, p3y :INTEGER; x, y, y1, y2, y3, lastX, nextX, thisX, thisY, corTop, corBot, vertexRadius, tmpReal :REAL; pt, begPt, endPt, lastPt, nextPt :POINT; totVol, correction, abundance, avgDen, decOfTotVol, totRad :REAL; correctPolyArray :DYNARRAY OF POINT; correctPolyCnt :INTEGER; pMinCor :INTEGER; t :REAL; bestPt1, bestPt2, bestPt3 :POINT; upperCore, lowerRadiative, upperRadiative, lowerConvective :REAL; scLower, scUpper, scDense, scTotal, previousScore :REAL; parallel, intOnLines :BOOLEAN; iterations :INTEGER; testMaxX, testMaxY, testMinX, testMinY :INTEGER;   {$Include ..\..\..\Common\Utilities_All.px} {$Include ..\..\..\Common\Utilities_Objects.px}   PROCEDURE CubicBezier(precision :INTEGER); VAR cnt :INTEGER; t, t2, t3 :REAL; BEGIN FOR cnt := 1 TO precision DO BEGIN t := cnt / precision; bezier[cnt] := {p0 * ((1 - t) ^ 3) p0 is always 0,0 for us, so skip it} + (p1 * ((1 - t) ^ 2) * 3 * (t )) + (p2 * ((1 - t) ) * 3 * (t ^ 2)) + (p3 * (t ^ 3)); END; END;   FUNCTION RadiusOfSphere(volume :REAL) :REAL; BEGIN RadiusOfSphere := (volume / ((4 * PI) / 3)) ^ (1/3); END;   PROCEDURE Calculate(p1x, p1y, p2x, p2y, p3y :INTEGER); BEGIN iterations := iterations + 1; IF (p1x + p1y + p2x + p2y + p3y) > 0 THEN BEGIN p1.x := p1x; p1.y := p1y; p2.x := p2x; p2.y := p2y; p3.y := p3y; CubicBezier(bezPrec); END;   {Add the correction factors to the abundances, and calculate the volumes and the radii.} totVol := 0; FOR cnt := 0 TO lastIndex DO BEGIN abundance := 10 ^ (rows[cnt].use + rows[cnt].cor); abundance := abundance * 100000000; rows[cnt].vol := abundance / rows[cnt].apv; totVol := totVol + rows[cnt].vol; totRad := RadiusOfSphere(totVol); rows[cnt].rad := totRad; END;   {totRad is now the total radius, so now we can convert the rads to decimals of the total. And knowing the total volume, we can calculate the average density.} upperCore := 0; lowerRadiative := 0; upperRadiative := 0; lowerConvective := 0; avgDen := 0; FOR cnt := 0 TO lastIndex DO BEGIN rows[cnt].rad := rows[cnt].rad / totRad; IF (rows[cnt].rad < .27 - .05) THEN upperCore := rows[cnt].den; IF (rows[cnt].rad < .70 - .05) THEN upperRadiative := rows[cnt].den; IF ((rows[cnt].rad > .27 + .05) AND (lowerRadiative = 0)) THEN lowerRadiative := rows[cnt].den; IF ((rows[cnt].rad > .70 + .05) AND (lowerConvective = 0)) THEN lowerConvective := rows[cnt].den; decOfTotVol := rows[cnt].vol / totVol; avgDen := avgDen + (decOfTotVol * rows[cnt].den); END; scDense := Abs(1408 - avgDen) / -1000; scLower := (upperCore - lowerRadiative) / 2000; scUpper := (upperRadiative - lowerConvective) / 2000; scTotal := scDense + scLower + scUpper; IF (scTotal > previousScore) THEN BEGIN bestPt1 := p1; bestPt2 := p2; bestPt3 := p3; previousScore := scTotal; END; END;   PROCEDURE PlaceText(text :STRING; x, y, rot :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); END;   BEGIN IF GetCustomObjectInfo(objName, objHand, recHand, wallHand) THEN BEGIN SetObjectVariableBoolean(objHand, 800, TRUE); IF pZoom THEN zoomFactor := 5 ELSE zoomFactor := 1;   {Load the abundance data.} Open(dataFile); FOR cnt := 0 TO lastIndex DO BEGIN ReadLn( rows[cnt].sym, rows[cnt].nam, rows[cnt].den, rows[cnt].apv, rows[cnt].set1, rows[cnt].set2, rows[cnt].set3); END; Close(dataFile); FOR cnt := 0 TO lastIndex DO BEGIN IF pSet = '1' THEN rows[cnt].use := rows[cnt].set1 ELSE IF pSet = '2' THEN rows[cnt].use := rows[cnt].set2 ELSE IF pSet = '3' THEN rows[cnt].use := rows[cnt].set3; END;   IF (debug) THEN ReWrite(debugOutFile); IF (debug) THEN Close(debugOutFile); OpenPoly;   p0.x := 0; p0.y := 0; p1.x := pControlPoint01X; p1.y := pControlPoint01Y; p2.x := pControlPoint02X; p2.y := pControlPoint02Y; p3.x := xMax; p3.y := pControlPoint03Y; SetRField(objHand, objName, 'ControlPoint03X', Concat(xMax)); IF ((pReCalc) and (false)) THEN BEGIN iterations := 0; previousScore := -100; testMinX := -10; testMinY := -50; testMaxX := 200; testMaxY := 50; FOR p1x := testMinX TO testMaxX DO BEGIN FOR p1y := testMinY TO testMaxY DO BEGIN FOR p2x := testMinX TO testMaxX DO BEGIN FOR p2y := testMinY TO testMaxY DO BEGIN FOR p3y := testMinY TO testMaxY DO BEGIN Calculate(p1x, p1y, p2x, p2y, p3y); END; END; END; END; END; p1 := bestPt1; p2 := bestPt2; p3 := bestPt3; Message('iterations: ', iterations); SetRField(objHand, objName, 'ControlPoint01X', Concat(p1.x)); SetRField(objHand, objName, 'ControlPoint01Y', Concat(p1.y)); SetRField(objHand, objName, 'ControlPoint02X', Concat(p2.x)); SetRField(objHand, objName, 'ControlPoint02Y', Concat(p2.y)); SetRField(objHand, objName, 'ControlPoint03Y', Concat(p3.y)); SetRField(objHand, objName, 'ReCalc', 'False'); END; CubicBezier(bezPrec); IF pShow_Bezier THEN BEGIN BeginPoly; AddPoint(0, 0); FOR cnt := 1 TO bezPrec DO BEGIN AddPoint(bezier[cnt].x, bezier[cnt].y); END; EndPoly; SetFPat(LNewObj, 0); SetPenFore(LNewObj, 65535, 0, 0); END;   {Get the vertices of the correction polyline.} pathHand := GetCustomObjectPath(objHand); ReshapablePath(objHand, pathHand); cnt := GetVertNum(pathHand); GetPolylineVertex(pathHand, cnt, x, y, vertexType, vertexRadius); SetPolylineVertex(pathHand, cnt, xMax + .001, y, vertexType, vertexRadius, TRUE); pathHand := CopyPoly(pathHand); SetFPat(pathHand, 0); SetVertexVisibility(pathHand, GetVertNum(pathHand) - 1, FALSE); correctPoly := ConvertToPolygon(pathHand, 32); correctPolyCnt := GetVertNum(correctPoly); ALLOCATE correctPolyArray[1..correctPolyCnt]; FOR cnt := 1 TO correctPolyCnt DO BEGIN GetPolyPt(correctPoly, cnt, pt.x, pt.y); correctPolyArray[cnt] := pt; END; DelObj(correctPoly); IF NOT pZoom THEN BEGIN FOR cnt := 1 to GetVertNum(pathHand) DO BEGIN GetPolylineVertex(pathHand, cnt, x, y, vertexType, vertexRadius); SetPolylineVertex(pathHand, cnt, x, -y / 5, vertexType, vertexRadius, TRUE); END; SetLS(pathHand, -1); END;   {Find y values along correction polyline.} begPt.y := 100; endPt.y := -100; FOR cnt1 := 0 TO lastIndex DO BEGIN begPt.x := cnt1 * 2; endPt.x := cnt1 * 2; correction := 0; FOR cnt2 := 1 to correctPolyCnt - 1 DO BEGIN lastPt := correctPolyArray[cnt2]; nextPt := correctPolyArray[cnt2+1]; LineLineIntersection(begPt, endPt, lastPt, nextPt, parallel, intOnLines, pt); IF ((NOT parallel) & (intOnLines)) THEN BEGIN correction := correction + (pt.y / 5); cnt2 := correctPolyCnt; {quit loop} END; END; rows[cnt1].cor := correction; END;   Calculate(0, 0, 0, 0, 0);   SetRField(objHand, objName, 'Density', Concat(avgDen)); SetRField(objHand, objName, 'scDense', Concat(scDense)); SetRField(objHand, objName, 'scLower', Concat(scLower)); SetRField(objHand, objName, 'scUpper', Concat(scUpper)); SetRField(objHand, objName, 'scTotal', Concat(scTotal));   {***********} {Graphics...} {***********} pMinCor := pMaxCor - 12; corTop := pMaxCor * zoomFactor; corBot := pMinCor * zoomFactor;   {Abundances.} y1 := corTop + 23; IF NOT pZoom THEN y1 := y1 - 25; TextSize(5); FOR cnt := 0 TO lastIndex DO BEGIN {element values and names} nextPt.x := cnt * 2; nextPt.y := rows[cnt].use + y1; Arc(nextPt.x - .2, nextPt.y + .2, nextPt.x + .2, nextPt.y - .2, 0, 360); IF cnt > 0 THEN LineTo(nextPt.x, nextPt.y) ELSE MoveTo(nextPt.x, nextPt.y); PlaceText(rows[cnt].nam, nextPt.x + .4, y1 - 3, -90, 'left', 'center'); END; Rect(0, y1 - 2, xMax, y1 + 12); SetFPat(LNewObj, 0); FOR cnt := -2 TO 12 DO BEGIN {vertical tick marks} IF cnt MOD 2 = 0 THEN BEGIN MoveTo(0, y1 + cnt); LineTo( - 1, y1 + cnt); MoveTo(xMax, y1 + cnt); LineTo(xMax + 1, y1 + cnt); PlaceText(Concat(cnt), xMax + 3.5, y1 + cnt, 0, 'right', 'center'); PlaceText(Concat(cnt), - 1.5, y1 + cnt, 0, 'right', 'center'); END; END;   {Corrections.} IF pZoom THEN BEGIN Rect(0, corBot, xMax, corTop); SetFPat(LNewObj, 0); FOR cnt := pMinCor TO pMaxCor DO BEGIN {vertical tick marks} IF cnt MOD 2 = 0 THEN BEGIN MoveTo(0, cnt * zoomFactor); LineTo( - 1, cnt * zoomFactor); MoveTo(xMax, cnt * zoomFactor); LineTo(xMax + 1, cnt * zoomFactor); PlaceText(Concat(cnt), xMax + 3.5, cnt * zoomFactor, 0, 'right', 'center'); PlaceText(Concat(cnt), - 1.5, cnt * zoomFactor, 0, 'right', 'center'); END; END; END;   {y axis labels} TextSize(8); y3 := corBot - 50; FOR cnt := 0 TO 20 DO BEGIN MoveTo(0, y3 + cnt); LineTo( - 1, y3 + cnt); MoveTo(xMax, y3 + cnt); LineTo(xMax + 1, y3 + cnt); IF cnt MOD 2 = 0 THEN BEGIN IF cnt MOD 4 = 0 THEN BEGIN MoveTo(0, y3 + cnt); LineTo( - 2, y3 + cnt); MoveTo(xMax, y3 + cnt); LineTo(xMax + 2, y3 + cnt); PlaceText(Concat(cnt), xMax + 6.5, y3 + cnt + .1, 0, 'right', 'center'); PlaceText(Concat(cnt), - 2.5, y3 + cnt + .1, 0, 'right', 'center'); END; END; END;   {density curve} MoveTo(0, (rows[0].den / 1000) + y3); FOR cnt := 1 TO lastIndex DO BEGIN LineTo(xMax * rows[cnt].rad, (rows[cnt-1].den / 1000) + y3); LineTo(xMax * rows[cnt].rad, (rows[cnt ].den / 1000) + y3); END;   {Export results to a text file.} WriteLn('0.0000000000', Chr(9), rows[0].den); FOR cnt := 0 TO lastIndex DO BEGIN WriteLn(rows[cnt].rad, Chr(9), rows[cnt].den); IF (cnt < lastIndex) THEN WriteLn(rows[cnt].rad + 0.0000000001, Chr(9), rows[cnt + 1].den); END;   {primary elements} TextSize(9); MoveTo(0, y3 + 22); LineTo(thisX, y3 + 22 + 1); lastX := 0; FOR cnt := 0 TO lastIndex DO BEGIN thisX := xMax * rows[cnt].rad; thisY := y3 + 22; MoveTo(thisX, thisY); LineTo(thisX, thisY + 1); IF (thisX - lastX > 6) THEN BEGIN PlaceText(rows[cnt].sym, (thisX + lastX) / 2, thisY + .25, 0, 'center', 'bot'); END; lastX := thisX; END;   FOR cnt := 0 TO 10 DO BEGIN MoveTo((cnt / 10) * xMax, y3 - 2); LineTo((cnt / 10) * xMax, y3 - 3); END;   IF pZoom THEN BEGIN TextSize(8); PlaceText(Concat('resulting average density: ', Round(avgDen), ' kg/m '), xMax * .500, corBot - 1, 0, 'center', 'top'); TextSize(6); PlaceText('3', (xMax * .500) + 27.7, corBot - 1, 0, 'center', 'top'); END;   Rect(0, y3 - 2, xMax, y3 + 22); SetFPat(LNewObj, 0); MoveTo(xMax * .27, y3 - 2); LineTo(xMax * .27, y3 + 22); SetLS(LNewObj, -1); MoveTo(xMax * .70, y3 - 2); LineTo(xMax * .70, y3 + 22); SetLS(LNewObj, -1);   TextSize(12); IF pSet = '1' THEN str := ' in Photosphere, from Anders & Grevesse (1989)' ELSE IF pSet = '2' THEN str := ' in Solar Wind, from Manuel, Kamat, & Mozina (2006)' ELSE IF pSet = '3' THEN str := ' in B2FH Values, from Manuel, Kamat, & Mozina (2006)'; PlaceText(Concat('Log(Abundance)', str), xMax * .500, y1 + 13, 0, 'center', 'bot'); PlaceText('Density per Solar Radius', xMax * .500, y3 + 27, 0, 'center', 'bot'); PlaceText('Core', xMax * .135, y3 - 3, 0, 'center', 'top'); PlaceText('Radiative Zone', xMax * .485, y3 - 3, 0, 'center', 'top'); PlaceText('Convective Zone', xMax * .850, y3 - 3, 0, 'center', 'top');   TextSize(10); PlaceText('Mg/m', 0, 0, 0, 'center', 'bot'); SetTextOrientation(LNewObj, -7, y3 + 9, 90, false); TextSize( 8); PlaceText( '3', 0, 0, 0, 'center', 'top'); SetTextOrientation(LNewObj, -12, y3 + 15, 90, false);   IF NOT pZoom THEN BEGIN tmpReal := 17.5; thisX := (xMax / 2) - (192.04 / 2); thisY := y1 + 25; Rect(thisX, thisY, thisX + 192.04, thisY - 123.46 + 25); SetFPat(LNewObj, 0); END; END; END; RUN(Main);