;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;/home/nakajima/heliog/cdaw10/gyrokspect4_3.pro
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;goto,jump1
;Input parameters ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
ex = 3.7
A = 2.7d+32*5.0 ;A*E^-ex (E in MeV)  ;10
B = 121.4        ;Gauss
theta1 = 40  ;45.0    ;degree
area = 27.0*14.0  ;arcsec^2 
Ds = 19.0  ;depth in arcsec
;energy 10^(0.1*j-2) MeV;j-0~40
 ;0-0.01,7-0.05,10-0.1,26-4.0,27-5.0,30-10,35-31.6,37-50
j1 = 7           ;lowest energy
j2 = 29 ;highest energy
;
;alpha = 3.0      ;Razin par,1.5fB/fp  ;20
Thot = 2.2e+7
Nhot = 3.0e+10
Tcold = 1.0e+4
Ncold = 1.5e+5  ;neglect cold temp plasma
Namb = Nhot+Ncold    ;ambient density
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;fpB = 1.5/alpha    ;=fp/fB
;fp = fpB*fB
;Namb= (fp/8.98e+3)^2 
fB = 2.8e+6*B                ;30
fp = 8.98e+3*sqrt(Namb)
fpB = fp/fB
alpha = 1.5/fpB
D = Ds*7.25e+7    ;depth in cm
V = area*Ds*3.81e+23  ;volume in cm^3
;non-thermal number density
NT = A/V   ;effective  ;;;;;;;;;;;;;;;;;;;;
theta = theta1*0.01745     ;rad
omega = area/4.2545e+10    ;sterad
;                            ;40
;Pitch-angle distribution ;;;;;;;;;;;;;;;;;;;;;;
g = 1.0/(4.0*3.1416)
;
hnumb = fltarr(301)
freq = fltarr(301)
;
;ordinary mode
Gfo = dblarr(301,41)
Hfo = dblarr(301,41)
Go1 = dblarr(301,40)        ;50
Go = dblarr(301)
Ho1 = dblarr(301,40)
Ho = dblarr(301)
jo1 = dblarr(301)   ;emissivity/NT*B
ko1 = dblarr(301)   ;abs coefficient/(NT/B)
jo = dblarr(301)    ;emissivity
ko = dblarr(301)    ;emissivity
;
Gfo(*,*)=0
Hfo(*,*)=0                 ;60
Go1(*,*)=0
Go(*)=0
Ho1(*,*)=0
Ho(*)=0
jo1(*)=0
ko1(*)=0
jo(*)=0
ko(*)=0
;
;extraordinary mode        ;70
Gfx = dblarr(301,41)
Hfx = dblarr(301,41)
Gx1 = dblarr(301,40)
Gx = dblarr(301)
Hx1 = dblarr(301,40)
Hx = dblarr(301)
jx1 = dblarr(301)   ;emissivity/NT*B
kx1 = dblarr(301)   ;abs coefficient/(NT/B)
jx = dblarr(301)    ;emissivity
kx = dblarr(301)    ;abs coefficient    ;80
;
Gfx(*,*)=0
Hfx(*,*)=0
Gx1(*,*)=0
Gx(*)=0
Hx1(*,*)=0
Hx(*)=0
jx1(*)=0
kx1(*)=0
jx(*)=0                     ;90
kx(*)=0
;
de = fltarr(40)
flo = dblarr(301)
flx = dblarr(301)
S = dblarr(301)
poldg = fltarr(301)
;
eff = dblarr(301)  ;f-f emissivity
kff = fltarr(301)  ;f-f absorption coefficient   ;100
floff = dblarr(301)
flxff = dblarr(301)
Sff = dblarr(301)  ;include f-f absorption
Sff2 = dblarr(301)  ;include f-f absorption
poldgff = fltarr(301)
;
de(*)=0
flo(*)=0
flx(*)=0
S(*)=0                  ;110
Sff(*)=0
Sff2(*)=0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
for i=0,300 do begin 
hn = 10.0^(0.01*i)   ;min=1,max=1000
hnumb(i) = hn
freq(i) = hn*fB
;
;refraction index ;;;;;;;;;;;;;;;;;;;;;;;;;;;
nord = (1.0+2.0*(fpB)^2*(fpB^2 - hn^2)/$    ;120
 ((hn^4*(sin(theta))^4+4.0*hn^2*(fpB^2-hn^2)^2*(cos(theta))^2)^0.5$
  -2.0*hn^2*(fpB^2-hn^2)-hn^2*(sin(theta))^2))^0.5
next = (1.0+2.0*(fpB)^2*(fpB^2 - hn^2)/$
 (-(hn^4*(sin(theta))^4+4.0*hn^2*(fpB^2-hn^2)^2*(cos(theta))^2)^0.5$
  -2.0*hn^2*(fpB^2-hn^2)-hn^2*(sin(theta))^2))^0.5
;
;polarization coefficient ;;;;;;;;;;;;;;;;;;;;;;;;
atheto = -2.0*hn*(fpB^2 - hn^2)*cos(theta)/(-hn^2*(sin(theta))^2$
  +(hn^4*(sin(theta))^4+4.0*hn^2*(fpB^2 - hn^2)^2*(cos(theta))^2)^0.5)
athetx = -2.0*hn*(fpB^2 - hn^2)*cos(theta)/(-hn^2*(sin(theta))^2$  ;130
  -(hn^4*(sin(theta))^4+4.0*hn^2*(fpB^2 - hn^2)^2*(cos(theta))^2)^0.5)
;
;flux density and abs. coefficient of free-free radiation ;;;;;;;;;
eff(i) = 3.1d-37*freq(i)*Thot*omega*(1.0-exp(9.78e-3*D*Nhot^2/$
  freq(i)^2/Thot^1.5*(24.5+alog(Thot)-alog(freq(i)))))
kff(i) = 9.78e-3*Ncold^2/(freq(i)^2)/(Tcold^1.5)$
  *(17.7+alog(Tcold^1.5)-alog(freq(i)))   ;*0.54?
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
for j =j1, j2-1 do begin              ;140
el = 10^(0.1*j-2)
eu = 10^(0.1*(j+1)-2)
em = 10^(0.1*(j+0.5)-2)
de(j) = eu-el
gam = em/0.511 +1.0
beta = sqrt(1.0-gam^(-2))
;
;for ordinary mode
myuo = nord*beta*cos(theta)
so = gam*hn*(1.0-nord*beta*myuo*cos(theta))   ;150
xo =nord*beta*sqrt(1.0-myuo^2)*sin(theta)/$
  (1.0-nord*beta*myuo*cos(theta))
zo = xo*exp(sqrt(1.0-xo^2))/(1.0+sqrt(1.0-xo^2))
ao = ((1.0-xo^2)^1.5+0.503297/so)^0.1666667
bo = ((1.0-xo^2)^1.5+1.193000/so)^0.1666667*(1.0-0.2/so^0.6666667)
fo = ao*bo
epo = 1.0/sqrt(1.0-nord^2*beta^2)
tauo = epo*beta*nord*sin(theta)
;
;for extraordinary mode              ;160
myux = next*beta*cos(theta)
sx = gam*hn*(1.0-next*beta*myux*cos(theta))
xx =next*beta*sqrt(1.0-myux^2)*sin(theta)/$
  (1.0-next*beta*myux*cos(theta))
zx = xx*exp(sqrt(1.0-xx^2))/(1.0+sqrt(1.0-xx^2))
ax = ((1.0-xx^2)^1.5+0.503297/sx)^0.1666667
bx = ((1.0-xx^2)^1.5+1.193000/sx)^0.1666667*(1.0-0.2/sx^0.6666667)
fx = ax*bx
epx = 1.0/sqrt(1.0-next^2*beta^2)
taux = epx*beta*next*sin(theta)      ;170
;
kapp = (beta^2*gam^2*(ex*gam+2.0*(gam-1.0))+(gam-1.0)) $
  /(beta^2*gam^3*(gam-1.0))
u = em^(-ex)
;
Gfo(i,j) = sqrt(3.1416*hn)*(-fo*(1.0+tauo^2)+atheto*cos(theta))^2 $
 /(nord*(sin(theta))^2*(1.0+atheto^2)*ao^2)*u*g*zo^(2.0*so) $
 /(sqrt(gam)*epo^3*(1.0+tauo^2)^0.25)
;
Hfo(i,j) = sqrt(3.1416*hn)*(-fo*(1.0+tauo^2)+atheto*cos(theta))^2 $ ;180
 /(nord*(sin(theta))^2*(1.0+atheto^2)*ao^2)*u*g*zo^(2.0*so) $
 /(sqrt(gam)*epo^3*(1.0+tauo^2)^0.25)/hn^2*kapp/nord
;
Gfx(i,j) =sqrt (3.1416*hn)*(-fx*(1.0+taux^2)+athetx*cos(theta))^2 $
 /(next*(sin(theta))^2*(1.0+athetx^2)*ax^2)*u*g*zx^(2.0*sx) $
 /(sqrt(gam)*epx^3*(1.0+taux^2)^0.25)
;
Hfx(i,j) = sqrt(3.1416*hn)*(-fx*(1.0+taux^2)+athetx*cos(theta))^2 $
 /(next*(sin(theta))^2*(1.0+athetx^2)*ax^2)*u*g*zx^(2.0*sx) $
 /(sqrt(gam)*epx^3*(1.0+taux^2)^0.25)/hn^2*kapp/next         ;190
;
endfor
;
;integration ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
Go1(i,j1)=de(j1)*Gfo(i,j1)   ;dblarr(40)
Ho1(i,j1)=de(j1)*Hfo(i,j1)
Gx1(i,j1)=de(j1)*Gfx(i,j1)
Hx1(i,j1)=de(j1)*Hfx(i,j1)
;                               ;200
for k = j1+1,j2-1 do begin
;
Go1(i,k) = Go1(i,k-1)+de(k)*Gfo(i,k)
Ho1(i,k) = Ho1(i,k-1)+de(k)*Hfo(i,k)
Gx1(i,k) = Gx1(i,k-1)+de(k)*Gfx(i,k)
Hx1(i,k) = Hx1(i,k-1)+de(k)*Hfx(i,k)
;
endfor
;
;G and H function          ;210
Go(i) = Go1(i,j2-1)
Ho(i) = Ho1(i,j2-1)
Gx(i) = Gx1(i,j2-1)
Hx(i) = Hx1(i,j2-1)
;
jo1(i) = (1.354e-22)*Go(i)    ;NT*B: neglected
ko1(i) = (1.896e-8)*Ho(i)  ;NT/B: neglected
jx1(i) = (1.354e-22)*Gx(i)    ;NT*B: neglected
kx1(i) = (1.896e-8)*Hx(i)  ;NT/B: neglected
;                          ;220
;Emissivity & abs coefficient  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;
jo(i) = jo1(i)*NT*B
ko(i) = ko1(i)*NT/B
jx(i) = jx1(i)*NT*B 
kx(i) = kx1(i)*NT/B
;
;Flux density;;;; eff(i),flux density of f-f emision, can be added
flo(i) = jo(i)/ko(i)*(1.0-exp(-ko(i)*D))*omega   ;erg/(s,Hz,cm^2)
flx(i) = jx(i)/kx(i)*(1.0-exp(-kx(i)*D))*omega   ;erg/(s,Hz,cm^2)
S(i) = (flo(i)+flx(i))*1.0e+19   ;eff(i);sfu      ;230
poldg(i)= (flx(i)-flo(i))/(flx(i)+flo(i))
;
;Flux density including f-f absorption  ;;;;;;;;;;;;;;;;;;;
;
floff(i) = jo(i)/(ko(i)+kff(i))*(1.0-exp(-(ko(i)+kff(i))*D))*omega   
flxff(i) = jx(i)/(kx(i)+kff(i))*(1.0-exp(-(kx(i)+kff(i))*D))*omega
;;;;eff(i),flux density of f-f emision, can be added;
Sff(i) = (floff(i)+flxff(i))*1.0e+19   ;eff(i)
poldgff(i) = (flxff(i)-floff(i))/(floff(i)+flxff(i))
;                                                  ;240
Sff2(i) = (jo(i)+jx(i))/(ko(i)+kx(i)+kff(i))*$
 (1.0-exp(-(ko(i)+kx(i)+kff(i))*D))*omega*1.0e+19  ;wrong equation?
;
endfor
;
;print,kff(198)*D
;
;window,0,xs=512,ys=512
;plot_oo,hnumb,ko1,xr=[1.0,1000.0]
;plot_oo,hnumb,kx1,xr=[1.0,1000.]       ;250
;plot_oo,freq,S,xr=[1.0e+9,1.0e+11]
;
;jump1:
;goto,jump3
;
freq2 = [2.e+9,3.75e+9,9.4e+9,1.7e+10,3.5e+10,8.0e+10]
flux = [720.0,1370.0,2770.0,1730.0,810.0,138.0]
Nfreq = n_elements(freq2)
errors = replicate(0.0, Nfreq)
errors[Nfreq-1] = 0.4  ;0.8                  ;260
;
;jump3: 
;goto,jump4
set_plot,'ps'
device,/portrait,xs=18.0,ys=9.0,xoffset=2.0,yoffset=3.0
device,file='gyrokspect4_3.eps',/encaps
;plot_oo,hnumb,jx1,xr=[1.0,100.]
;oplot,hnumb,jo1
;oplot,hnumb,kx1,linestyle=2
;oplot,hnumb,ko1,linestyle=2            ;270
;plot_oo,freq,S,xr=[1.0e+9,1.0e+12],yr=[1.0,1.0e+5]
plot_oo,freq,S,xr=[1.0e+9,1.0e+11],yr=[10.0,1.0e+4],/xstyle,$
 thick=3.0,linestyle=0,xtitle='Frequency (GHz)' $
 ,ytitle='Flux Density (SFU)' $
 ,title='Radio Spectrum obtained from HXR One (Peak 1 of Event 5)'
;oplot,freq,Sff,linestyle=0,thick=3.0
;;plot_oi,freq,poldg,xr=[1.0e+9,1.0e+11],yr=[-1.0,1.0],/xstyle
;
oploterr,freq2, flux, errors*flux, psym = 4,thick=3,symsize =1.0  ;0.6
;
device,/close                           ;280
set_plot,'x'
;
;jump4:
;
spind=-3.2257*alog10(S(201)/S(170))  ;34.784/17.036, B=121.4G
;spind=-2.7778*alog10(S(237)/S(201))   ;79.685/34.7846, B=121.4
;spind=-3.227*alog10(S(231)/S(200))    ;34.7/17, B=60.7
print,S(170),S(201),spind
;
end