print, 'Program for plotting FFT spectra of NCO-A.' print, 'Run this program after "ncoapow.pro"' print, ' ' print, 'Output Medium: ' print, ' ' print, ' output.ps (0)' print, ' output.eps (1)' print, ' Tektronix-Screen (2)' print, ' X-Screen (3)' print, ' Windows (4)' print, ' terminate program (5)' read, medium case medium of 0: begin set_plot, 'PS' device, filename='output.ps', $ bits = 8, /color, $ /landscape, $ xoffset=2.5, yoffset=26., $ xsize=24., ysize=16. ; /portrait, $ ; xoffset=2., yoffset=3., $ ; xsize=16., ysize=24. end 1: begin set_plot, 'PS' device, filename='output.eps', $ bits = 8, /color, $ /encapsulated, $ xsize=24., ysize=16. end 2: begin set_plot, 'TEK' end 3: begin set_plot, 'X' end 4: begin set_plot, 'WIN' end else: begin print, 'Programmende' stop end endcase black = 0 white = 1 brightness = findgen(10)/10. one = replicate(1,10) RED = [black, white, one, brightness, 0.5+brightness/2., one, one, brightness, one] GREEN = [black, white, brightness, one, 0.5+brightness/2., 0.5+brightness/2., brightness, one, one] BLUE = [black, white, brightness, brightness, one, brightness, one, one, brightness] tvlct, 255 * RED, 255 * GREEN, 255 * BLUE red = 2 & green = 12 & blue = 22 & orange = 32 & magenta = 42 & cyan = 52 & yellow = 62 xyouts, 0.5, 0.5, /normal, '!5' dummy = ' ' xmin = 0.0 xmax = 0.0 print, 'frequency range:' print, 'Enter minimum frequency for display (-1 for 0-4 Hz):' read, xmin if xmin ne -1 then begin print, 'Enter maximum frequency for display:' read, xmax endif else begin xmin = 0.0 xmax = 4.0 endelse ;ymin = 0.0 ;ymax = 0.0 ;print, 'power range:' ;print, 'Enter minimum power for display (-1 for 10^-1 to 10^6 Hz):' ;read, ymin ;if ymin ne -1 then begin ; print, 'Enter maximum power for display:' ; read, ymax ;endif else begin ; ymin = 1.e-1 ; ymax = 1.e6 ;endelse ymin = 1.e-1 ymax = 1.e6 print, 'Read data? (y/n) read, dummy if dummy eq 'y' then begin fname = ' ' ptit = ' ' read, fname, prompt = 'Enter file name for NCO-A power: ' read, ptit, prompt = 'Enter name and date of test (plot title): ' read, dummy, prompt = 'File ' + fname + ' covers how many minutes? ' ; ptit = ptit + ': NCO-A FFT over 150 min' ptit = ptit + ': NCO-A FFT over ' + dummy + ' min' openr, 1, fname ln = long(0.0) while not eof(1) do begin readf, 1, dummy ln = ln + 1 endwhile close, 1 powdat = dblarr(2, ln) openr, 1, fname readf, 1, powdat close, 1 index_fosc = where(powdat(0,*) gt 0.365 and powdat(0,*) lt 0.369) fosc = powdat(0, index_fosc(0) + where(powdat(1,index_fosc) eq max(powdat(1,index_fosc)))) print, powdat(*,index_fosc) print, 'oscillation frequency = ' + string(fosc, format = '(f6.4)') + ' Hz' print, 'hit any key to accept oscillation frequency or "m" to modify' dummy = get_kbrd(1) if dummy eq 'm' then $ read, fosc, prompt = 'Enter modifyed value for oscillation frequency: ' index_fo3 = where(powdat(0,*) gt 0.002 and powdat(0,*) lt 0.005) fo3 = powdat(0, index_fo3(0) + where(powdat(1,index_fo3) eq max(powdat(1,index_fo3)))) print, powdat(*,index_fo3) repeat begin print, 'beat frequency = ' + string(fo3, format = '(f6.4)') + ' Hz' fo2 = (fosc + fo3) / 2. print, 'resulting frequency fo2 = (fosc + fo3) / 2.0' print, 'fo2 = ' + string(fo2, format = '(f6.4)') + ' Hz' print, 'hit any key to accept beat frequency or "m" to modify' dummy = get_kbrd(1) if dummy eq 'm' then $ read, fo3, prompt = 'Enter modifyed value for beat frequency: ' endrep until dummy ne 'm' endif !p.charsize = 1.5 !p.charthick = 2.0 plot_io, powdat(0,*), powdat(1,*), $ color = black, /nodata, $ position = [0.1, 0.0, 1.0, 0.9], $ ; title = fname, $ ; title = ptit, $ ;see below xtitle = 'FREQUENCY [Hz]', $ ytitle = 'POWER [Hz!E2!N/Hz]', $ xrange = [xmin, xmax], $ xstyle = 9, $ xticklen = -0.02, $ yrange = [ymin, ymax], $ ystyle = 1, $ xtickname = replicate('',30) step = 1. i_max = 4. ;j_max = 8. j_max = 10. k_max = 6. for i = 0., i_max do $ ;for i = -i_max, i_max do $ for j = -j_max, j_max do $ for k = -k_max, k_max do begin xpos = i + j * fosc + k * fo2 if xpos(0) ge xmin or xpos(0) le xmax then $ oplot, [xpos, xpos], [ymin, ymax/10.*10.^step], $ color = ((abs(i) eq 0.0) * red + (abs(i) eq 1.0) * magenta + (abs(i) eq 2.0) * green + (abs(i) eq 3.0) * cyan + (abs(i) eq 4.0) * blue) + (abs(j) + abs(k)) * 9. / (j_max + k_max), $ thick = 4.0 - 3.8 * abs(j) / j_max, $ linestyle = (abs(k) eq 0.0) * 0 + (abs(k) eq 1.0) * 5 + (abs(k) eq 2.0) * 2 + (abs(k) eq 3.0) * 3 + (abs(k) eq 4.0) * 4 + (abs(k) ge 5.0) * 1 step = step - 1. / (i_max + 1.) / (2. * j_max + 1.) / (2. * k_max + 1.) ; let lines become shorter in each step to see overplotting endfor oplot, powdat(0,*), powdat(1,*), color = black, thick = 1.0 axis, xmin, ymax/10., /xaxis, xstyle = 1, $ color = black, $ xrange = [xmin / fosc, xmax / fosc], $ xtickname = replicate('',30), $ ; xtitle = 'Multiples of f!D1!N = ' + string(fosc, format = '(f6.4)') + ' Hz' xtitle = '!AMultiples of !18f!N1!X' axis, /xaxis, xstyle = 1, $ color = black, $ xtickname = replicate('',30), $ xrange = [(4.0 - xmin) / fosc, (4.0 - xmax) / fosc] xyouts, 0.5, 1.0, /normal, $ ptit, alignment = 0.5, $ charsize = 2.0 poly1x = [xmin+(xmax-xmin)/40., xmin+(xmax-xmin)/5., xmin+(xmax-xmin)/5., xmin+(xmax-xmin)/40., xmin+(xmax-xmin)/40.] poly1y = [ymax/100.*5., ymax/100.*5., ymax/100., ymax/100., ymax/100*5.] polyfill, poly1x, poly1y, $ color = white oplot, poly1x, poly1y, thick = 0.5, color = black xyouts, xmin+(xmax-xmin)/9., ymax/100.*10.^0.45, alignment = 0.5, $ '!18f!D1!N!X = ' + string(fosc, format='(d6.4)') + ' Hz', $ charsize = 1.0, charthick = 1.0 xyouts, xmin+(xmax-xmin)/9., ymax/100.*10.^0.15, alignment = 0.5, $ '!18f!D2!N!X = ' + string(fo2, format='(d6.4)') + ' Hz', $ charsize = 1.0, charthick = 1.0 ;xyouts, xmin+(xmax-xmin)/9., ymax/100.*10.^0.15, alignment = 0.5, $ ; '!18f!D3!N!X = ' + string(fo3, format='(d6.4)') + ' Hz', $ ; charsize = 1.0, charthick = 1.0 !p.charsize = 1.0 !p.charthick = 1.0 if (medium eq 0) or (medium eq 1) then begin device, /close spawn, 'mv -f output.ps ncoapow_' + string(xmin, format = '(d5.3)') + '-' + string(xmax, format = '(d5.3)') + '.ps' endif end