print, 'calculates FFT-spectrum of NCO-A' dummy = ' ' ;---- t_start and t_end determination --------------------------- openr, 1,'carriera.dat' datc = dblarr(2) readf, 1, datc t_start = datc(0) while not eof(1) do readf, 1, datc close, 1 t_end = datc(0) ;----- Abfragen -------------------------------------------------- fname = 'ncoa.dat' print, 'Time scale:' print, ' ' print, ' t_0 = 0 (1)' print, ' Start of test = 0 (2)' print, ' scale with given time (3)' read, tscale case tscale of 1: begin print, 'Enter moment of t_0 as recorded from HK3 (in hours)' read, t_sub end 2: begin t_sub = t_start end else: begin t_sub = 0.0 end endcase t_start = t_start - t_sub t_end = t_end - t_sub print, 'Which time range for FFT ?' print, 'Enter minimum time in minutes (-1 for full range)' read, mintime if mintime ne -1 then begin print, 'Enter maximum time in minutes' read, maxtime endif else begin mintime = 0.0 maxtime = (t_end - t_start) * 60.0 endelse res_flag = 0 print, 'Use NCO-A residuals ? (y/n)' print, '(a 10-sec-smooth will be subtracted from NCO-A data)' read, dummy if dummy eq 'y' then res_flag = 1 ;---- Daten lesen --------------------------------------------------- print, 'read data ? (y/n)' read, dummy if dummy eq 'y' then begin openr, 1, fname ln = long(0.0) while not eof(1) do begin readf, 1, dummy ln = ln + 1 endwhile close, 1 dat = dblarr(2,ln) dat_res = dblarr(2,ln) openr, 1, fname readf, 1, dat close, 1 dat(0,*) = (dat(0,*) - t_sub) * 60.0 ; Minuten dat_res(0,*) = dat(0,*) dat_res(1,*) = dat(1,*) - smooth(dat(1,*), 81) endif index_start = long(-1.0) index_end = long(-1.0) i = long(0.0) while (index_start lt 0) do begin if dat(0,i) ge mintime then index_start = i i = i + 1 endwhile while (index_end lt 0) do begin if dat(0,i) ge maxtime then index_end = i i = i + 1 endwhile nop = index_end - index_start + 1 ;--------FFT-transform---------------------------------- t_fft = dindgen(nop) * 8.0 / (nop - 1) if res_flag eq 0 then $ f_fft = fft(dat(1,index_start:index_end), 1) $ else $ f_fft = fft(dat_res(1,index_start:index_end), 1) power = 2.0 * abs(f_fft) * abs(f_fft) / nop / 8.0 openw, 1, 'ncoapow.dat' i = long(0.0) while t_fft(i) le 4.0 do begin printf, 1, format = '(d20.10, e20.10)', t_fft(i), power(i) i = i + 1 endwhile close, 1 ;-------Bestimmung des Osc.-Amplitude------------------ index_osc = where(t_fft gt 0.35 and t_fft lt 0.38) maxpow = max(power(index_osc), index_max) print, ' ' print, 'Maximum Power: ' + string(maxpow) + ' Hz2/Hz' print, 'at f = ' + string(t_fft(index_max+index_osc(0))) + ' Hz' print, 'Resolution: ' + string(t_fft(1)-t_fft(0)) + ' Hz' end