Source code :: main

[Return]
  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
151
;---------------------------
;IDL MAKEFILE
;---------------------------
;
;---------------------------
;Central file of used IDL analysis
;---------------------------

pro main                                                       ;has to correspond to the name specified in call command in makefile

;---------------------------
;DAS PARAMETERS INPUT
;---------------------------

dt_slow = 0.5D-5
dt_fast = 1d-6
polarity = [1d,1d,-1d,-1d,1d,-1d,1d,1d,1d,-1d,-1d,1d,1d,1d,0d,-1d]  ;[-]
normalisation = [1d,1d,1d,1d,1d,2d,2d,2d,2d,2d,2d,2d,2d,2d,0d,1d]   ;[-]
uloop_cal = 5.5d  ;[-]
bt_cal = 70.42d   ;[T/Wb]

;---------------------------
;INITIALISATION 
;---------------------------

path_fast = '../../../../DAS/0311NIturbo.ON/'+['NIturbo_01','NIturbo_02','NIturbo_03','NIturbo_04','NIturbo_05','NIturbo_06','NIturbo_07','NIturbo_08','NIturbo_09','NIturbo_10','NIturbo_11','NIturbo_12','NIturbo_13','NIturbo_14','NIturbo_15','NIturbo_16']           ;location of central data with respect to folder containing this program
path_slow = '../../0411ShotHomepage.ONN/usbscopes'   ;location of basic parameters data
read = size(file_reading(4,path_slow),/dimensions)
N_slow = read[1]
read = size(file_reading(2,path_fast[0]), /dimensions)
N_fast = read[1]
shot_no = long(file_reading(1, '../../../../ShotNo'))
time_slow = dindgen(N_slow)*dt_slow*1d3 ;[ms]
time_fast = dindgen(N_fast)*dt_fast*1d3 ;[ms]
theta = 360d/16d*(dindgen(15)-1d) ;[deg] poloidal angle

plasmadetect = long(file_reading(1,'../../0411ShotHomepage.ONN/Plasma'))
if plasmadetect eq 1L then begin
plasmastart = file_reading(1,'../../0411ShotHomepage.ONN/PlasmaStart')
plasmaend = file_reading(1,'../../0411ShotHomepage.ONN/PlasmaEnd')
plasmastart = plasmastart*1d3 ;[mus]
plasmaend = plasmaend*1d3     ;[mus]
time_low = plasmastart[0]
time_up = plasmaend[0]
endif else begin
time_low = 8d3
time_up = 16d3
endelse
N_fft = 400d                  ;[mus] length of fft window
fft_jmp_rat = 0.10d           ;[-] ratio of time_len by what the fft windows are moved forward
N_fft_time = long((time_up-time_low)/N_fft/fft_jmp_rat+1d)
a_w = 2d                      ;[mus] coeff. in STFT for Gaussian windowing
N_vid_time = 200d
vid_jmp_rat = 0.1d
N_vid = long((time_up-time_low)/N_vid_time/vid_jmp_rat + 1)

;---------------------------
;BASIC PARAMETERS 
;---------------------------

;FILE READING
data_slow = file_reading(4,path_slow)

;CALCULATION
offset_1 = total(data_slow[0,0:499])/5d2
offset_2 = total(data_slow[1,0:499])/5d2
offset_3 = total(data_slow[2,0:499])/5d2
data_slow[0,*] = uloop_cal*(data_slow[0,*]-offset_1)
data_slow[1,*] = bt_cal*total(data_slow[1,*]-offset_2,/cumulative)*dt_slow
data_slow[2,*] = 1d-3*i_plasma(data_slow[0,*],data_slow[2,*]-offset_3,dt_slow)
q_est = 90.3125d*data_slow[1,*]/data_slow[2,*]

;PLOTTING
make_plot = plot_multi(time_slow,[data_slow[0:2,*],q_est],['Uloop [V]','Bt [T]','Ip [kA]','q [-]'],'Basic parameters of shot no. '+strtrim(shot_no,2),'../outputs/plot_'+strtrim(shot_no,2)+'_basic',[-1d,-5d-2,-2d-1,0d],[15d,0.5d,6d,8d])

;---------------------------
;MHD ANALYSIS 
;---------------------------

;FILE READING AND DATA PREPARATION
data_fast = dblarr(16,N_fast)
for i=0,15 do begin
data_temp = file_reading(2,path_fast[i])
data_fast[i,*] = data_temp[1,*]
endfor
for i=0,15d do begin
data_fast[i,*] = polarity[i]*normalisation[i]*data_fast[i,*]
endfor
for i=0,15 do begin
data_fast[i,*] = smooth(data_fast[i,*],5)
endfor

;FFT ANALYSIS
FFT_arr = dblarr(N_fft_time, N_fft/2, 16)
fft_time = time_low + dindgen(N_fft_time)/double(N_fft_time)*(time_up-time_low)  ;[mus]
fft_freq = dindgen(N_fft/2)/double(N_fft/2)/dt_fast/2d
data_fft = data_fast
w_t = exp(-0.5d*(a_w*(dindgen(N_fft)-0.5d*double(N_fft))/(0.5d*double(N_fft)))^2)  ;windowing function (Gaussian)
for i=0,15 do begin
  for j=0, N_fft_time-1 do begin
    data_temp = abs(fft(w_t*data_fft[i,time_low+double(j)*fft_jmp_rat*N_fft:time_low+(double(j)*fft_jmp_rat+1d)*N_fft-1], /double))
    FFT_arr[j,*,i] = smooth(data_temp[0:N_fft/2-1],3)
  endfor
endfor
for i=0,8 do begin
  make_plot = plot_fft(fft_time, 1d-3*fft_freq[0:N_fft/7-1], fft_arr[*,0:N_fft/7-1,i],'Time [ms]','f [kHz]','Frequencies at coil no. '+strtrim(i+1,2)+' of shot no. '+strtrim(shot_no,2), '../outputs/plot_'+strtrim(shot_no,2)+'_fft_0'+strtrim(i+1,2))
endfor
for i=9,15 do begin
  make_plot = plot_fft(fft_time, 1d-3*fft_freq[0:N_fft/7-1], fft_arr[*,0:N_fft/7-1,i],'Time [ms]','f [kHz]','Frequencies at coil no. '+strtrim(i+1,2)+' of shot no. '+strtrim(shot_no,2), '../outputs/plot_'+strtrim(shot_no,2)+'_fft_'+strtrim(i+1,2))
endfor

;VIDEO PLOTS
data = dblarr(15,N_fast)              ;this is all due to coil 15 being off
data[0,*] = data_fast[15,*]
data[1:14,*] = data_fast[0:13,*]
if long(N_vid) le 99L then begin
	for i=0,9 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_000'+strtrim(i,2),0.5d)
	endfor
	for i=10,N_vid-1 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_00'+strtrim(i,2),0.5d)
	endfor
endif
if (long(N_vid) gt 99L) and (long(N_vid) le 999L) then begin
	for i=0,9 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_000'+strtrim(i,2),0.5d)
	endfor
	for i=10,99 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_00'+strtrim(i,2),0.5d)
	endfor
	for i=100,N_vid-1 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_0'+strtrim(i,2),0.5d)
	endfor
endif
if (long(N_vid) gt 999L) then begin
	for i=0,9 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_000'+strtrim(i,2),0.5d)
	endfor
	for i=10,99 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_00'+strtrim(i,2),0.5d)
	endfor
	for i=100,999 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_0'+strtrim(i,2),0.5d)
	endfor
	for i=1000,N_vid-1 do begin
		make_plot = plot_profile(time_fast[time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1], theta, transpose(data[*,time_low+double(i)*vid_jmp_rat*N_vid_time:time_low+(double(i)*vid_jmp_rat+1d)*N_vid_time-1]),'Time [ms]','Theta [deg]','Perturbations of shot no. '+strtrim(shot_no,2), '../outputs/video/plot_pert_'+strtrim(i,2),0.5d)
	endfor
endif


end

Navigation