ecr_dispersion/ecr_dispersion.pro at master · dlg0/ecr_dispersion · GitHub

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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
;+
;	Calculate the kPer or kPar for waves
;	in the electron cyclotron range based
;	on Brambilla.
;
;	David L Green, ORNL, 28-Jul-2010
;	greendl1@ornl.gov
;-


pro ecr_dispersion, $
	kPar = kPar, $
	kPer = kPer, $
	SI = SI

	if (not keyword_set(kPar)) $
			and (not keyword_set(kPer)) then begin

		print, 'Input error:'
		print, 'Either kPer or kPar must be set.'
		print, 'For kPer=0 or kPar=0 just use 1e-6 or the like ;)'

	endif


	; parameters
	; ----------

	f	= 28.0d9	; Hz


	if not keyword_set ( SI ) then begin

		; Stupid Brambilla units (cgs)
		; ----------------------------

		print, ''
		print, '*** Using CGS units'
		print, '*** Ensure inputs are appropriate'
		print, ''

		ne_	= 1.0d12	; cm^-3
		b0	= 1.1d4		; Gauss = 1d4*Tesla

		print, 'ne: ', ne_, ' [cm^-3]'
		print, 'B: ', b0, ' [Gauss]'
		print, ''

		e 	= 4.8032d-10
		Z	= -1d0
		q	= e * Z
		me	= 9.109d-28
		c	= 2.99792d10

		wce	=	q * b0 / ( me * c )
		wpe	= sqrt ( 4 * !pi * ne_ * q^2 / me )

		units = ' [cm^-1]'

	endif else begin

		; SI units
		; --------

		print, ''
		print, '*** Using SI units'
		print, '*** Ensure inputs are appropriate'
		print, ''

		ne_	= 1.0d18	; m^-3
		b0	= 1.1 		; T

		print, 'ne: ', ne_, ' [m^-3]'
		print, 'B: ', b0, ' [T]'
		print, ''

		e 	= 1.602d-19
		Z	= -1d0
		q	= e * Z
		me	= 9.109d-31
		e0 	= 8.85d-12
		c	= 2.99792d8

		wce = q * b0 / me
		wpe	= sqrt ( ne_ * q^2 / ( me * e0 ) )

		units = ' [m^-1]'

	endelse

	w	= 2 * !pi * f


	; From pg 203-207 Brambilla
	; -------------------------

	X 	= wpe^2 / w^2
	u 	= wce / w

	if keyword_set(kPar) then begin

		nPar = kPar * c / w

		del1	= u^2 * ( 2 - nPar^2 )^2 + 4 * nPar^2 * ( 1-X )

		if del1 lt 0 then begin
			print, 'Error: del1 < 0'
			stop
		endif

		nPerSq_O	= 1 - nPar^2 - X $
				- X*u/2 * ( u*(1+nPar^2)-sqrt(del1) ) / ( 1 - X - u^2 )
		nPerSq_X	= 1 - nPar^2 - X $
				- X*u/2 * ( u*(1+nPar^2)+sqrt(del1) ) / ( 1 - X - u^2 )

		print, 'kPer^2_O: +/-', nPerSq_O*w^2/c^2
		print, 'kPer^2_X: +/-', nPerSq_X*w^2/c^2

		print, 'kPer_O: +/-', sqrt(nPerSq_O*w^2/c^2), units
		print, 'kPer_X: +/-', sqrt(nPerSq_X*w^2/c^2), units

	endif

	if keyword_set(kPer) then begin

		nPer = kPer * c / w

		del2	= ( 2 * ( 1-X ) - nPer^2 )^2 + nPer^4 * ( u^2 - 1 )

		if del2 lt 0 then begin
			print, 'Error: del2 < 0'
			stop
		endif


		nParSq_L	= 1.0 / ( 2 * ( 1-u^2 ) ) $
			   	* ( 1-X-u^2+X*u * ( nPer^2*u + sqrt(del2) ) / (1-X) )
	   	nParSq_R	= 1.0 / ( 2 * ( 1-u^2 ) ) $
				* ( 1-X-u^2+X*u * ( nPer^2*u - sqrt(del2) ) / (1-X) )

		print, 'nPar^2_L: +/-', nParSq_L*w^2/c^2
		print, 'nPar^2_R: +/-', nParSq_R*w^2/c^2

		print, 'kPar_L: +/-', sqrt(nParSq_L*w^2/c^2), units
		print, 'kPar_R: +/-', sqrt(nParSq_R*w^2/c^2), units

	endif


end