;;;;;;;;;;;;;;;;;; Read Data Files ;;;;;;;;;;;;;;;;;;;;;;; Pro read_dat, shotn, file, channel, data, Mistake ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; !Quiet = 1 @ida3_setup filn = strcompress(string(shotn), /remove_all) Mistake = 0 shot = '00'+strmid(filn,0,2)+'.'+strmid(filn,2,2) ; filn = strmid(shot, 2, 2)+strmid(shot, 5, 2) filp = strcompress(file+shot, /remove_all) ;filename = strcompress('/net/fuslsa/data/MAST_RES/'+filn+'/'+filp, /remove_all) filename = strcompress('/net//fuslsa/data/MAST_Data/'+filn+'/LATEST/'+filp, /remove_all) filenum = Ida_open(filename, ida_read) If (filenum eq 0) then begin print, 'File: ', filename, ' Not Found' print, 'shot:', shot print, 'file:', file print, 'chan:', channel Mistake = 1 Goto, jumpMistake Endif Item_err = Ida_find(filenum, itemname=channel) If (item_err eq 0) then begin print, 'Data Item: ', channel, ' Not Found' Mistake = 1 Goto, jumpMistake Endif err = Ida_tvec(item_err, t) err = Ida_get_data(item_err, y, ida_d4+ida_real) tmax = Max(t, Imax) data = fltarr(2, Imax+1) data(1,*) = y data(0,*) = t err = ida_free(item_err) err = ida_close(filenum) jumpmistake : !quiet = 0 End ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Calibrate, dat ;;;;;;;;;;;;;;;;;;;Calibration from shot #7798 C1= [0.726572, 1.31234, 1.26832, 0.911739, 1.15886, 1.17279,$ 0.704816, 1.11000, 1.11520, 1.041983, 1.17517, 0.958079,$ 0.973466, 1.07078, 0.837652, 1.15311, 1.15231, 0.815988,$ 1.08168, 0.898918, 1.01480, 1.06571, 0.883064, 0.769014,$ 1.04758, 1.23981, 1.40572, 1.10865, 0.915037, 1.12535,$ 1.21061, 1.28835] C2= [0.48, 0.59, 0.574, 0.712, 0.516, 0.964, $ 1.012, 0.89, 1.096, 0.72, 0.719, 0.938, $ 0.953, 1.18, 3.22, 12.35, 23.11, 9.21, $ 3.07, 1.185, 0.905, 0.810, 1.166, 1.2, $ 1.495, 0.838, 1.043, 1.199, 2.08, 2.91, $ 2.805, 3.057] For i=0,31 Do Begin dat(*,i)=dat(*,i)*C1(i)*10 dat(*,32+i)=dat(*,32+i)*C2(i)*10 EndFor Return, dat End ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro ebw_emiss,nshot,t_slice, Sm, Nf read_dat,nshot,'xne', 'XNE_EBW1', EB1, error if error lt 0 then begin print,"error reading XNE_EBW1,error=",error return end read_dat,nshot,'xne', 'XNE_EBW2', EB2, error if error lt 0 then begin print,"error reading XNE_EBW2,error=",error return end read_dat,nshot,'xne', 'XNE_EBW3', EB3, error if error lt 0 then begin print,"error reading XNE_EBW3,error=",error return end read_dat,nshot,'xne', 'XNE_EBW4', EB4, error if error lt 0 then begin print,"error reading XNE_EBW4,error=",error return end read_dat,nshot,'xne', 'XNE_EBW5', EB5, error if error lt 0 then begin print,"error reading XNE_EBW5,error=",error return end read_dat,nshot,'xne', 'XNE_EBW6', EB6, error if error lt 0 then begin print,"error reading XNE_EBW6,error=",error return end read_dat,nshot,'xne', 'XNE_EBW7', EB7, error if error lt 0 then begin print,"error reading XNE_EBW7,error=",error return end read_dat,nshot,'xim', 'xim_da/bo6*', Dal, error if error lt 0 then begin print,"error reading xim_da/bo6*,error=",error return end read_dat,nshot,'amc', 'amc_plasma current', Ip, error if error lt 0 then begin print,"error reading amc_plasma current,error=",error return end read_dat,nshot,'ane', 'ane_density', N_e, error if error lt 0 then begin print,"error reading ane_density,error=",error return end f = [15.88, 16.32, 17.04, 17.94, 18.90, 19.86, 20.84, 21.76, $ 22.68, 23.60, 24.38, 25.10, 25.76, 26.38, 26.94, 27.48, $ 25.52, 26.16, 27.00, 27.98, 29.22, 30.50, 31.74, 32.96, $ 34.18, 35.38, 36.46, 37.48, 38.38, 39.22, 39.94, 40.60] sf = fltarr(64) For i = 0,31 Do sf(2*i) = f(i) For i = 0,14 Do Begin sf(2*i+1) = (f(i)+f(i+1))*0.5 sf(2*i+33) =(f(i+16)+f(i+17))*0.5 ; sf(2*i+65) =(f(i+32)+f(i+33))*0.5 EndFor sf(31)=27.75 sf(63)=40.9 ; sf(95)=61.40 ;;;;;;;;;; Define regular time time = EB7(0,*) n = floor(N_elements(time)/32.) timestep = (Max(time)-Min(time))/n tr = Min(time)+timestep*findgen(n) ;;;;;;;;;;;;;;;;;;;;;;;;;Convert Log Scale;;;;;;;;;;;;;;;;;; EB1(1,*) = 1.97*exp(EB4(1,*)/0.34-9.)+0.03 ;EB3(1,*) = 2*exp(EB5(1,*)/0.34-9.) EB2(1,*) = 2.73*exp(EB6(1,*)/0.34-9.)+0.014 ;;;;;;;;;;;;;;;;;;;;;;;;;Decode signals;;;;;;;;;;;;;;;;;;;;; lim = -.05 datt = fltarr(n, 64) For i=0,31 Do Begin si = where(EB7(1,*) GT lim AND EB7(1,*) LT lim+0.153/2) dat = EB1(1,si) t = EB1(0,si) dt = interpol(dat,t,tr) dn = total(dt(0:20))/21. datt(*,i) = (dt - dn) dat = EB2(1,si) t = EB2(0,si) dt = interpol(dat,t,tr) dn = total(dt(0:20))/21. datt(*,i+32) = (dt - dn) ; dat = EB3(1,si) ; t = EB3(0,si) ; dt = interpol(dat,t,st) ; dn = total(dt(0:10))/11. ; datt(*,i+64) = Abs(dt - dn) lim = lim+.153/2 EndFor datt = Calibrate(datt) ;;;;;;;;;;;;;;;;;;;;;;;;;;;Bad frequency filtering;;;;;;;;;;;;;;;;;; Bf = fltarr(64) For i=0,63 Do Bf(i) = stddev(datt(*,i)) print, Bf goodi = where((Bf GT Nf(0)) and (Bf LT Nf(1))) datt=datt(*,goodi) sf = sf(goodi) print, 'Number of frequencies=', n_elements(sf) ;;;;;;;;;;;;;;;;;;;;;;;;;;;Sorting by frequency and Smoothing ;;;;;; find = sort(sf) f = sf(find) datt = datt(*,find) If Sm GE 1 Then datt=smooth(datt,Sm) ;;;;;;;;;;;;;;;;;;;;;;;;; Making regular matrix ;;;;;;;;;;;;;;;; fn = N_elements(f) stepf = (Max(f)-Min(f))/(fn-1) startf = Min(f) fr = startf + stepf*findgen(fn) datr = fltarr(n, fn) For j = 0,n-1 Do datr(j,*) = interpol(datt(j,*),f,fr) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; tn = N_elements(t_slice) ts = replicate(1.,tn,64) For i=0,tn-1 Do ts(i,*)=t_slice(i)*ts(i,*) tin = intarr(tn) Device, Decomposed=0 ;!P.Color = 25 !P.Multi=[0,1,2,0,0] Loadct,27 Contour, datr, tr, fr, /Follow, NLevels=29, /Fill, $ Charsize=1, background=1, $ ZRange=[0, 0.8], $ XRange=[t_slice(0), t_slice(tn-1)], $ Xstyle=1, Ystyle=1 , $ Title = 'EBW Spectrum #'+string(nshot),$ Xtitle = 'Time, ms', $ Ytitle = 'Frequency, GHz' ;;;;;;;;;;;;;;;;;;;;; Color Setting ;;;;;;;;;;;;;;;;;;;;;;;;; red = [0,1,1,0,0,1,1,0] green=[0,1,0,1,0,0,1,1] blue =[0,1,0,0,1,1,0,1] TVLCT,255*red,255*green,255*blue ;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plot,[t_slice(0), t_slice(tn-1)],[Min(f),Max(f)],/Nodata,/Noerase, $ ; Xstyle=1, Ystyle=1, background=0, $ ; Title = 'EBW Spectrum #'+string(nshot),$ ; Xtitle = 'Time, ms', $ ; Ytitle = 'Frequency, GHz', $ ; Charsize=2, color=0 ; For i=0,tn-1 Do Oplot, ts(i,*),f, color=i+2 For i=0,tn-1 Do Begin mn = Min(Abs(tr-t_slice(i)),ind) tin(i)=ind Endfor Plot, [15,41], [0,2.0],/Nodata,$ ; Title='EBW spectrum',$ Thick = 300,$ Xtitle='Frequency, GHz',$ Ytitle='EBW signal, a.u.' For i=0,tn-1 Do Oplot, fr, datr(tin(i),*), color=i+2 ; Plot, Dal(0,*), Dal(1,*), color=0, background=1,$ ; YRange=[0, 1.0], $ ; XRange=[t_slice(0), t_slice(tn-1)], $ ; Xstyle=1, Ystyle=1, $ ; Charsize=2,$ ; Title='D_alpha, Red:*10^20, Blue:Ip,MA',$ ; Xtitle='Time, ms',$ ; Ytitle='Da,Ne,Ip, a.u.' ; Oplot, Dal(0,*), Dal(1,*), color=0 ; Oplot, Ip(0,*), Abs(Ip(1,*))*1e-3, color=4 ; Oplot, N_e(0,*),N_e(1,*)*2.1e-21, color=2 !P.Multi=0 ;;;;;;;;;;;;;;;;;;;;; Writing Data ;;;;;;;;;;;;;;;;;;;;;;;;; Openw,unit,'ebw_emiss.dat',/get_lun ; Printf,unit,'Frequency, GHz', ' Time=',strmid(strcompress(string(t_slice(0)),/remove_all),0,5),$ ; ' Time=',strmid(strcompress(string(t_slice(1)),/remove_all),0,5), ' shot#',nshot For i=0,fn-1 Do Begin For j=0,n-1 Do Printf,unit,tr(j), fr(i), datr(j,i) EndFor Close,unit Free_lun,unit Openw,unit,'ebw_spec.dat',/get_lun For i=0,fn-1 Do Begin Printf,unit,fr(i), datr(tin(1),i), datr(tin(2),i),datr(tin(3),i) EndFor Close,unit Free_lun,unit END