earthorbaxis.pro
3.17 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
;==================================================
; DD
; ORBITAL OBJECTS
; mexorbaxis.pro
; V.2.0
; Modifications:
; 19 May 2006: V.1.0, Fedorov
; 09 Dec 2006: V.2.0, Fedorov, New approach for graphical environment
;---------------------------------------------------
; Implementation of all types of axis for orbital plot
; Options:
; Key ToPlot
; /CYL - cylindrical XR 0
; /XY 1
; /XZ 2
; /YZ 3
;====================================================
pro earthorbaxis, ToPlot
common GraphC, graph, GraphN
;------------ Set graphics Colors ------------------------
if (Graph[GraphN].Lx.range[0] eq Graph[GraphN].Lx.range[1]) then begin
Graph[GraphN].Lx.range[0] = -20;
Graph[GraphN].Lx.range[1] = 10;
endif
if (Graph[GraphN].Ly.range[0] eq Graph[GraphN].Ly.range[1]) then begin
Graph[GraphN].Ly.range[0] = ToPlot eq 0 ? 0: -20;
Graph[GraphN].Ly.range[1] = 10;
endif
!x = Graph[GraphN].Lx
!y = Graph[GraphN].Ly
x = !x.range
y = !y.range
black = !DNC- 1
zero = fix(!DNC * 20.0/ 256.0 )
planet = fix(!DNC * 10.0/ 256.0)
bound = fix(!DNC * 1.0/ 256.0)
;------------------------------------------------------------
ThickPlanet = 2
ThickBound = 2
;------- main axis -----------------
case ToPlot of
0: begin
!x.title = 'X, R!DE'
!y.title = 'R, R!DE'
end
1: begin
!x.title = 'X, R!DE'
!y.title = 'Y, R!DE'
end
2: begin
!x.title = 'X, R!DE'
!y.title = 'Z, R!DE'
end
3: begin
!x.title = 'Y, R!DE'
!y.title = 'Z, R!DE'
end
endcase
;------ Plot main axis. Use common Graph system variables
plot,x,y,/nodata,/noerase,xstyle=1,ystyle=1,color=black,xminor=1, yminor=1, charsize=graph[GraphN].Lp.charsize
;------- plot Zero lines ---------------------
oplot,x,[0,0],color=zero, thick = 1
oplot,[0,0],y,color=zero, thick = 1
;------- plot Earth ---------------------
dal = 2.0 * !PI / 60.
if ToPlot EQ 0 then $
al = dal * findgen(31) $
else al = dal * findgen(61)
xx = cos(al)
yy = sin(al)
oplot,xx,yy,color=black, thick = ThickPlanet
if ToPlot EQ 0 then begin
;---- Plot boundaries for cylindrical plot
;----------------------- Bow Shock -----------------------------
RamPress = 1.5
theta = (135.*!pi/180.)*findgen(100)/100.
; rr = 2.04/(1.0+1.02*cos(theta))
xx = fltarr(100)
yy = fltarr(100)
epsilon = 0.81
L0 = 24.8
ramP = 1.8
L_BS = L0*(RamPress/ramP)^(-1./6.)
rr=L_BS/(1.+epsilon*COS(theta))
xx = rr *cos(theta)
yy = rr*sin(theta)
oplot,xx,yy,color=bound, thick = ThickBound
oplot,xx,-yy,color=bound, thick = ThickBound
;----------------------- MP Shue-----------------------------
Bz = -1.0
alpha=(.58-.007*Bz)*(1.+0.024*alog(RamPress))
r0=(10.22+1.29*tanh(0.184*(Bz+8.14)))*RamPress^(-1./6.6)
rr=r0*(2./(1.+COS(theta)))^alpha
xx = rr*cos(theta)
yy = rr*sin(theta)
oplot,xx,yy,color=bound, thick = ThickBound
oplot,xx,-yy,color=bound, thick = ThickBound
endif
return
end