;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;Program to calculate trapping electrons from injected
; electron flux (amp*E^delta), amp:triangular shape)
;/home/nakajima/heliog/trapsimul.pro 
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;goto,jump1
dr=500   ;burst duration in sec
m = 250  ; 1 sec step; usually dr/2?    ;10
dt = dr/(2.0*m)
amp = 100.0
damp = amp/m
delt = -3.7
namb10 = 1.0
p = 0.0
eng = fltarr(251)
;
;
ntrap = fltarr(501,251)   ;20
ntrap(*,*) = 0
;
;;Case1 for 0 <time< dr/2 
for j=1,250 do begin
  for k=0,250 do begin
    eng(k) = 20.0*k + 10.0
      for i=0,j do begin
;
  ntrap(j,k) = ntrap(j,k) + (1.0-p)*(damp*i)*(eng(k)^(delt)) $
    *dt*exp(-3.8*namb10*(j-i)*dt/eng(k))   ;30
;
      endfor
  endfor
endfor
;
;;Case2 for dr/2 <time <dr
for j=251,500 do begin
  for k=0,250 do begin
    eng(k) = 20.0*k + 10.0
;                           ;40
   for i=0,250 do begin
     ntrap(j,k) = ntrap(j,k) + (1.0-p)*(damp*i)*  $
     (eng(k)^(delt))*dt*exp(-3.8*namb10*(j-i)*dt/eng(k))
   endfor
   for i=251,j do begin
     ntrap(j,k) = ntrap(j,k) + (1.0-p)*(damp*(500-i))*  $
     (eng(k)^(delt))*dt*exp(-3.8*namb10*(j-i)*dt/eng(k))
   endfor
;
endfor           ;50
endfor
;
spectind=fltarr(501)
time=fltarr(501)
spectind(0)=0
for i=1,500 do begin
spectind(i)= -alog10(ntrap(i,50)/ntrap(i,40))/alog10(1010.0/810.0)
endfor
for j=0,500 do begin
time(j)=j   ;60
endfor
;
;jump1:
;goto,jump2
;
init_device=!d.name
set_plot,'ps'
device,file='trapsimul_1_2.eps',/encaps,/color,bits=8,/cmyk
loadct,12
;!p.multi=[0,1,2]       ;70
;
n=[350,450]
;[50,150,250,350,450,500]
plot_oo,eng(*),ntrap(1,*),xrange=[10,5000],yrange=[1.0e-11,1.0],$
 xstyle=1,ystyle=1
for i=0,1 do begin
oplot,eng(*),ntrap(n(i),*)  ; 40, 60, 80, 100
endfor
;
;jump2:          ;80
goto,jump3
;
init_device=!d.name
set_plot,'ps'
device,file='trapsimul_2.eps',/encaps,/color,bits=8,/cmyk;
;
plot,time,spectind,xrange=[0,500],yrange=[0,5]
;
jump3:
;
device,/close
set_plot, init_device
;              ;90
;
end