;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Program to calculate trapping electrons from injected
; electron flux (amp*E^delta), amp:triangular shape)
; delt:time-dependent (until 250 s)
;/home/nakajima/heliog/trapsimul2.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
delt0 = -3.7   ;*
delthrd = 1.0   ;*
namb10 = 1.0
p = 0.0
eng = fltarr(251)
;
;                         ;20
ntrap = fltarr(501,251)
ntrap(*,*) = 0
delt = fltarr(501)    ;*
;
;;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
;                         ;30
  delt(i) = delt0 + delthrd*i/250.0
;
  ntrap(j,k) = ntrap(j,k) + (1.0-p)*(damp*i)*(eng(k)^(delt(i))) $  ;*
    *dt*exp(-3.8*namb10*(j-i)*dt/eng(k))
;
      endfor
  endfor
endfor
;
;;Case2 for dr/2 <time <dr    ;40
for j=251,500 do begin
  for k=0,250 do begin
    eng(k) = 20.0*k + 10.0
;
   for i=0,250 do begin
;
     delt(i) = delt0 + delthrd*i/250.0
     ntrap(j,k) = ntrap(j,k) + (1.0-p)*(damp*i)*  $
     (eng(k)^(delt(i)))*dt*exp(-3.8*namb10*(j-i)*dt/eng(k))
;                             ;50
   endfor
   for i=251,j do begin
;
    delt(i) = delt0 + delthrd
     ntrap(j,k) = ntrap(j,k) + (1.0-p)*(damp*(500-i))*  $
     (eng(k)^(delt(i)))*dt*exp(-3.8*namb10*(j-i)*dt/eng(k))
;
   endfor
;
endfor           ;60
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   ;70
endfor
;
;jump1:
goto,jump2
;
init_device=!d.name
set_plot,'ps'
device,file='trapsimul2_1.eps',/encaps,/color,bits=8,/cmyk
loadct,12
;!p.multi=[0,1,2]       ;80
;
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:          ;90
;goto,jump3
;
init_device=!d.name
set_plot,'ps'
device,file='trapsimul2_2.eps',/encaps,/color,bits=8,/cmyk;
;
plot,time,spectind,xrange=[0,500],yrange=[0,5]
;
;jump3:
;               ;100
device,/close
set_plot, init_device
;
;
end