constellations.tcl

卫星仿真软件 卫星仿真软件 卫星仿真软件

#
######################################################
#
#  SaVi by Robert Thurman ([e-mail]) and
#          Patrick Worfolk ([e-mail]).
#
#  Copyright (c) 1997 by The Geometry Center.
#  This file is part of SaVi.  SaVi is free software;
#  you can redistribute it and/or modify it only under
#  the terms given in the file COPYRIGHT which you should
#  have received along with this file.  SaVi may be
#  obtained from:
#  http://savi.sourceforge.net/
#  http://www.geom.uiuc.edu/locate/SaVi
#
######################################################
#
# constellations
#
# Procedures to:
#
#   - generate a constellation
#   - read a two line element file
#
# TLE reading and constellation generation are incomplete, and
# hard to get right. Use at own risk.
# These need considerable expansion and work.
#
# $Id: constellations.tcl,v 1.4 2005/02/04 17:58:50 lloydwood Exp $


proc walker {T P F alt inc lan_offset phase_offset} {

    # dangerously confused; this attempted to generate a Walker delta
    # constellation, but mixed up Walker and Ballard notation.
    # Users might have expected a Walker star instead;
    # the term 'Walker constellation' is meaningless, since
    # Walker's papers cover a very wide range of geometries.
    # This procedure is retired, and now does nothing by design.
    # An altered version is below.

}


# Rough approximation of a Ballard (Walker delta) constellation.
# note use of harmonic factor AND use of phase_offset.
# First three parameters passed (NP, P, m) are in Ballard's notation.

#    T            total number of satellites
#    P            number of planes
#    F            interplane phasing is (360/T)*F - harmonic factor
#    alt          altitude of satellites
#                 - using radius, which is independent of the size
#                   of the central body, would be more rigorous.
#    inc          inclination of orbits to the equator.
#    lan_offset   longitude of ascending node offset
#    phase_offset offset for phase of satellites in planes
#		   - possibly redundant and might be misused, but as
#		     an initial offset for everything it's okay.
#
# By default, only one orbit in each orbital plane will be displayed
# since all the satellites move on the same orbit.
#
# We may want to pass (or compute?) a mask elevation angle to be set
# via upvar.
#

proc ballard_rosette {T P F alt inc lan_offset phase_offset} {
        global PI params

        # simulation constants set elsewhere in SaVi - see params.tcl
	set MU $params(Mu)
	set RADIUS_OF_EARTH $params(Radius)

	set NUM_PLANES $P
	set SATS_PER_PLANE [expr $T/$P]
	set INTERPLANE_SPACING [expr 360/$NUM_PLANES]

	# setup orbital elements
	set a [expr $alt+$RADIUS_OF_EARTH]
	set e 0.0
	set omega 0.0

	# compute period of orbit
	set T_per [expr 2 * $PI * pow($a,1.5) / sqrt($MU)]

	satellites GV_BEGIN
	for {set j 0} {$j < $NUM_PLANES} {incr j} {
		set Omega [expr $j * $INTERPLANE_SPACING + $lan_offset]
		set plane_offset [expr $T_per*($j*$F/$T + $phase_offset/360.0)]
		for {set i 0} {$i < $SATS_PER_PLANE} {incr i} {
			set T [expr $T_per * $i / $SATS_PER_PLANE + $plane_offset]
			set n [satellites LOAD $a $e $inc $Omega $omega $T]
			if {$i > 0} {
				satellites ORBIT_SET $n 0
			}
		}
	}
	satellites GV_END
}

#
# TLE reading can't yet handle satellites with different epochs,
# making it useless for more than one satellite. Much work is needed!
#
# Procedure to read a file in NORAD two line element (TLE) format
#
proc tle_file_input {filename} {
    global PI params

    set now [split [exec date -u {+%D %T %Y %j %H %M %S}] " "]
    set date [lindex $now 0]
    set time [lindex $now 1]
    puts stderr "For all these satellites, time 0 is $date $time GMT."
    set year [string trimleft [lindex $now 2] 0]
    set day [string trimleft [lindex $now 3] 0]
    set hour [string trimleft [lindex $now 4] 0]
    set minute [string trimleft [lindex $now 5] 0]
    set second [string trimleft [lindex $now 6] 0]
    set epoch_now [expr $day+($hour+($minute+$second/60.0)/60.0)/24.0]
    set MU $params(Mu)

    set f [open "$filename" r]
    set line1 0
    set line2 0
    savi GV_BEGIN
    while {[gets $f line] >= 0} {
	# remove white space
	set line [string trim $line]
	if {([string length $line] == 69) &&
	    ([string index $line 0] == 1)} {
	    # first elements line
	    set epoch_year [string range $line 18 19]
	    if {$epoch_year < 50} {
		incr epoch_year 2000
	    } {
		incr epoch_year 1900
	    }
	    set epoch [string range $line 20 31]
	    # compute delta t
	    set dt [time_difference $epoch_year $epoch $year $epoch_now ]
	    set line1 1
	} elseif {([string length $line] == 69) &&
		  ([string index $line 0] == 2) && ($line1 == 1)} {
	    # second elements line
	    set inc [string range $line 8 16]
	    set ascending_node [string range $line 17 24]
	    set ecc 0.[string range $line 26 32]
	    set arg_perigee [string range $line 34 41]
	    set mean_anomaly [string range $line 43 50]
	    set mean_motion [string range $line 52 62]
	    set line2 1
	} else {
	    # line with name or junk
	    set line1 0
	    set line2 0
	}

	# if we have everything then write out a satellite!
	if {($line1 == 1) && ($line2 == 1)} {
	    # write out orbital elements
	    set period [expr 86400/$mean_motion]
	    set a [expr pow($period*sqrt($MU)/2.0/$PI,2.0/3.0)]
	    set mean_anomaly [expr fmod($mean_anomaly + \
					    $dt*$mean_motion*360.0, 360.0)]
	    set T [expr $period*$mean_anomaly/360.0]
	    savi LOAD $a $ecc $inc $ascending_node $arg_perigee $T

	    set line1 0
	    set line2 0
	}
    }
    savi GV_END
}


#
# Compute the number of days between the first and second dates
# If second date is later, then positive.  The yi are years and
# the di are (fractional) days.
#
proc time_difference {y1 d1 y2 d2} {
    set dt [expr $d2-$d1]
    while {$y2 > $y1} {
	set dt [expr $dt+[days_in_year $y1]]
	incr y1
    }
    while {$y2 < $y1} {
	set dt [expr $dt-[days_in_year $y2]]
	incr y2
    }
    return $dt
}

#
# Compute the number of days in a given year
#
proc days_in_year {y} {
    # if not a fourth year return 365
    if { $y % 4 } { return 365 }
    # if a fourth year, but not a 100th year return 366
    if { $y % 100 } { return 366 }
    # if a 100th year, but not a 400th year return 365
    if { $y % 400} { return 365 }
    # it is a 400th year so return 366
    return 366
}