Attachment 'insertrb.pl'
Download 1 #!/usr/local/bin/perl
2 # insertrb.pl v0.2 2013 Sat 26
3 #
4 # This computes the two pairs of delta V burns to transfer a NORAD
5 # object to m288 orbit for rocket bodies
6 #
7 # ./insertrb.pl < {three line element file} > {output file}
8
9 use Astro::SpaceTrack;
10 use Astro::Coord::ECI;
11 use Astro::Coord::ECI::TLE;
12 use Astro::Coord::ECI::TLE::Set;
13 use Astro::Coord::ECI::Utils qw(deg2rad rad2deg);
14
15 my $r288 = 12789.0 ; # km
16
17 my $pi = 4.0*atan2(1.0,1.0) ;
18 my $mu = 3.986004418e5 ; # km3/s2 (no oblate effects)
19 my $v288 = sqrt($mu/$r288) ;
20 my $leo = 2000+6370 ; # low earth orbit semimajor radius
21
22 my @cislist = ("BLOCK","BREEZE","DNEPR","FREGAT","STRELA","SL-","SS-","SSN-");
23
24 my $cntall = 0 ;
25 my $cntdeb = 0 ;
26 my $cntrb = 0 ;
27 my $cntleo = 0 ;
28 my $cnthigh = 0 ;
29 my $cntcis = 0 ;
30 my $cntsl6 = 0 ;
31
32 # array for delta V histograms
33 my $rbhist = (0) x 30 ; # all rocket bodies
34 my @leohist = (0) x 30 ; # rocket bodies in LEO
35 my @highhist = (0) x 30 ; # rocket bodies above LEO
36 my $cishist = (0) x 30 ; # CIS/Soviet rockets
37 my @sl6hist = (0) x 30 ; # SL-6 above LEO (Molniya)
38
39 my $dvstep = 0.5 ; # array bin size 500m/s
40 my $dvmax = 11.0 ; # max delta V binned, 10 km/s
41 my $dvinmax = 22 ;
42
43 # printf("\n%7.0f=r288%7.3f=v288\n", $r288, $v288 );
44
45 open( F , ">rb.deltav");
46
47 foreach $el ( Astro::Coord::ECI::TLE->parse( <> ) ) {
48
49 my $nam = $el->get('name'); # name
50 my $bod = $el->body_type(); # body type
51 my $apo = $el->apogee(); # apogee km
52 my $per = $el->perigee(); # perigee km
53 my $sem = $el->semimajor(); # semimajor axis km
54 my $ecc = $el->get('eccentricity'); # eccentricity
55 my $inc = $el->get('inclination'); # inclination radians
56 my $apg = $el->get('argumentofperigee'); # arg of perigee radians
57 my $v01 = sqrt( $mu*(1.0-$ecc**2)/$sem ); # orbit velocity
58
59 $cntall += 1 ;
60 $cntdeb += 1 if ( $nam =~ / DEB/ ) ;
61
62 if( $nam =~ /R\/B/ ) {
63
64 # printf( "\n%7.3f=i%7.3f=ap%7.3f=v0", $inc, $apg, $v01 );
65 # printf( "\n%7.0f=sem%7.3f=ecc", $sem, $ecc );
66
67 # radius at equatorial crossings (two places)
68 my $r11 = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 + $ecc*cos($apg) ) ;
69 my $r12 = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 - $ecc*cos($apg) ) ;
70
71 # printf( "%7.0f=ra%7.0f=rp%7.0f=r11%7.0f=r12\n",$apo,$per,$r11,$r12 );
72
73 # tangential orbit velocity
74 my $vt11 = $v01 * ( 1.0 + $ecc*cos($apg) ) ;
75 my $vt12 = $v01 * ( 1.0 - $ecc*cos($apg) ) ;
76
77 # tangential orbit velocity in plane
78 my $vh11 = $vt11 * cos($inc);
79 my $vh12 = $vt12 * cos($inc);
80
81 # tangential orbit velocity out of plane
82 my $vv11 = $vt11 * sin($inc);
83 my $vv12 = $vt12 * sin($inc);
84
85 # radial orbit velocity
86 my $vr11 = $ecc * $v01 * sin($apg) ;
87 my $vr12 = -$ecc * $v01 * sin($apg) ;
88
89 # printf( "%7.3f=vt1%7.3f=vt2%7.3f=vr1%7.3f=vr2\n",$vt11,$vt12,$vr11,$vr12);
90
91 # transfer orbit semimajor axis
92 my $sem21 = 0.5 * ( $r11 + $r288 ) ;
93 my $sem22 = 0.5 * ( $r12 + $r288 ) ;
94
95 # transfer orbit eccentricity
96 my $ecc21 = ( $r11/$sem21 ) - 1.0 ;
97 my $ecc22 = ( $r12/$sem22 ) - 1.0 ;
98
99 # transfer orbit velocity
100 my $v021 = sqrt( $mu * ( 1.0 - $ecc21**2 ) / $sem21 );
101 my $v022 = sqrt( $mu * ( 1.0 - $ecc22**2 ) / $sem22 );
102
103 # transfer orbit r1 velocity
104 my $vh21 = ( 1.0 - $ecc21 ) * $v021 ;
105 my $vh22 = ( 1.0 - $ecc22 ) * $v022 ;
106
107 # transfer orbit r288 velocity
108 my $vh31 = ( 1.0 + $ecc21 ) * $v021 ;
109 my $vh32 = ( 1.0 + $ecc22 ) * $v022 ;
110
111 # deltaV to transfer orbit
112 my $dv11 = sqrt( $vr11**2 + $vv11**2 + ($vh21-$v288)**2 );
113 my $dv12 = sqrt( $vr12**2 + $vv12**2 + ($vh22-$v288)**2 );
114
115 # deltaV transfer orbit to m288
116 my $dv21 = abs( $vh21-$v288 ) ;
117 my $dv22 = abs( $vh22-$v288 ) ;
118
119 # sort the values, lowest first
120 if ( ( $dv12 < $dv11 ) && ( $dv22 < $dv21 ) ) {
121 my $tmp = $dv11 ; $dv11 = $dv12 ; $dv12 = $tmp ;
122 $tmp = $dv21 ; $dv21 = $dv22 ; $dv22 = $tmp ;
123 }
124
125 my $x = " "; # extra character
126 $x = "+" if($dv21 > $dv22) ;
127
128 my $dvt = $dv11 + $dv21 ;
129 my $dvindex = int ($dvt/$dvstep) ;
130
131 $rbhist[$dvindex]+=1; $cntrb+=1; $x="O";
132
133 foreach my $cis ( @cislist ) {
134 if( $nam =~ $cis )
135 { $cishist[$dvindex]+=1; $cntcis+=1; $x="S"; }
136 }
137
138 if( $nam =~ /SL-6/ )
139 { $sl6hist[$dvindex]+=1; $cntsl6+=1; $x="M"; }
140
141 if( $sem < $leo )
142 { $leohist[$dvindex]+=1; $cntleo+=1; }
143 else
144 { $highhist[$dvindex]+=1; $cnthigh+=1; }
145
146 printf F ("%7.3f%7.3f%8.3f%7.3f %s %s\n",
147 $dv11, $dv21, $dv12, $dv22, $x, $nam);
148 }
149 }
150 close F ;
151
152 my $cntncis = $cntrb - $cntcis ;
153
154 open( F , ">rb.stat");
155 printf F ("\ntotal=%5d", $cntall );
156 printf F ( " debris=%5d ", $cntdeb );
157 printf F ( " rocket bodies=%5d\n", $cntrb );
158 printf F ( "LEO R/B=%4d ", $cntleo );
159 printf F ( " High R/B=%4d\n", $cnthigh );
160 printf F ( "CIS/Soviet=%4d", $cntcis );
161 printf F ( " Soyuz/Molniya=%4d", $cntsl6 );
162 printf F ( " NonCIS=%4d\n", $cntncis );
163 close F ;
164
165 open( F, ">rb.dat" );
166 printf F ( "delV ALL HIGH CIS SL-6\n" );
167
168 for( my $i = 0 ; $i <= $dvinmax ; $i++ ) {
169 printf F ("%5.1f%6d%6d%6d%6d\n", $i*$dvstep,
170 $rbhist[$i], $highhist[$i], $cishist[$i], $sl6hist[$i] );
171 }
172 close F ;
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.