Scilab script for drawing stable orbits for RemoteTech constellations

RT_stable_orbit.sce

// Author: Dirtcrusher
// Last edited: 28/05/2015
//
// This work is licensed under the Creative Commons Attribution 4.0 International License.
// To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
// or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.


clear

// INPUT DATA
centralRadius = 600   // Radius of the central body
antennaRange  = 5000  // Range of the antenna
maxAmount     = 16    // Maximum number of satellites
maxAltitude   = 10000 // Maximum altitude of orbit
NbPoints      = 500   // Amount of points to draw
caption = "Kerbin - Communotron 32"

// INIT GRAPH
clf(0)
scf(0)

// CALCULATE DATA
altitudes = linspace(0, maxAltitude, NbPoints)
stableAltitudes = zeros(maxAmount, NbPoints)
minAmount = 3

for N=minAmount:maxAmount
    for i=1:NbPoints

        // Check if the centralBody blocks line of sight for connexion
        if altitudes(i) + centralRadius > centralRadius / cos(%pi / N) then
                
            // Calculate distance between 2 consecutives satellites
            distBetween = altitudes(i) * sqrt(2 - 2 * cos(2 * %pi / N))
            
            // Check if connexion between 2 consecutives satellites is possible
            if distBetween < antennaRange then
            
                // At this moment, the satellite constellation is stable
                stableAltitudes(N,i) = (altitudes(i) + centralRadius) * cos(%pi / N) + sqrt((antennaRange ^ 2) - 0.25 * (distBetween^2)) - centralRadius
            end
        end
    end
    
    // If the current number of satellites doesn't suffice for a stable connexion, don't draw it
    if max(stableAltitudes(N,:)) == 0 then
        minAmount = N + 1
    end
end

// PLOTTING
drawlater()

// Setup some nice drawing
fig = gcf()
fig.color_map = jetcolormap(maxAmount)
fig.figure_name = caption
xgrid()
xtitle(caption)

// Plot the actual curves
for N=minAmount:maxAmount
    
    // Only plot the points where stable orbit is not null
    imin = find(stableAltitudes(N,:)<>0, 1) // Index of first non-null datapoint
    imax = find(stableAltitudes(N,imin:NbPoints)==0, 1) - 1 // Index of last non-null datapoint
    
    if imax > 0 then
        // Add some text at the upper limit
        xstring(altitudes(imax), stableAltitudes(N,imax) / 3, string(N) + " satellites")
        
        // Draw curve with starting and ending limit
        plot2d([altitudes(imin), altitudes(imin:imax), altitudes(imax)], [0, stableAltitudes(N,imin:imax), 0], style=N)
    
    else
        // If the stable orbit goes on beyond the graph, we only need to draw
        // the starting limit
        plot2d([altitudes(imin), altitudes(imin:NbPoints)], [0, stableAltitudes(N,imin:NbPoints)], style=N)
    end
    // Increase thickness of drawings (need to do this here, for each graph, 
    // or else the background grid will follow the thickness)
    entity = gce()
    entity.children.thickness = 2
end

// label axis
xlabel("Orbit altitude (km)")
ylabel("Stable orbit altitude (km)")

// Actually draw everything
drawnow()

// If you want to save the figure to a file, uncomment the next line
//xs2png(fig, caption + ".png")