!! !! FtsPlatformingCorrection4.class !! to correct 30m FTS spectra for platforming (200kHz or 50kHz resolution) !! Baselines (1.order) are fitted to individual FTS units and subtracted. !! It is possible to define line windows to be ignored by the baseline !! fitting. See text marked "user input". !! !! written by Carsten Kramer, Manuel Gonzalez, Christof Buchbender !! based on ideas of Jerome Pety, 2-Feb-2012 !! Last update: v4.0 14-Jun-2012 !! !! Example application to a single spectrum: !! file in FTSOdp20120611.30m !! set scan 128 128 !! find; list !! get 28 !! set mode y total; set mode x total; set unit v f !! @ FtsPlatformingCorrection4 ! to define procedures !! @ subband ! to run correction procedure !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Debug option: 1 shows plots during correction and pauses ! 0 does not paus if .not.exist(debug) then def int debug /global endif let debug 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! List of Lines with center velocity and width !!! which are expected to appear in the !!! bandwidth of the Spectra that has to be !!! corrected. !!! Different Receivers can be treated at once !!! Lines that are not in the bandwidth of the !!! actual spectrum are ignored. !!! You can add as many lines as you want !!! Or leave it entirly empty !!! the latter defaults to 'set window 0 0'. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Start of user input: begin data gag_scratch:lines.dat -534 100 end data gag_scratch:lines.dat !! End of user input. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Procedure to get the channel for a specific Frequency ! ! Input &1 Frequency in GHz !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! begin procedure freqToChannel ! read the Input define real freqInput let freqInput &1 define real diffFreqCentralFreq define int channelDiff ! variable 'chan' stores the corresponding channel ! it has to be global for later use in prodecure win if .not.exist(chan) then define int chan /global endif ! calculate the difference between central frequency ! (stored in R%HEAD%SPE%RESTF) and the input frequency let diffFreqCentralFreq freqInput-R%HEAD%SPE%RESTF ! get the number of channels corresponding to ! that difference let channelDiff diffFreqCentralFreq/freq_step ! add the difference in channels to the ! central channel let chan R%HEAD%SPE%RCHAN+channelDiff exa chan end procedure freqToChannel !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 ! Get the Velocity for a specific Frequency ! ! Input &1 Frequency in GHz !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! begin procedure freqToVelocity ! read the Input define real freqInput let freqInput &1 define real diffFreqCentralFreq define int veloDiff ! variable 'velo' stores the corresponding velocity ! it has to be global for later use in prodecure win if .not.exist(velo) then define int velo /global endif ! calculate the difference between central frequency ! (stored in R%HEAD%SPE%RESTF) and the input frequency let diffFreqCentralFreq freqInput-R%HEAD%SPE%RESTF ! get the difference in velocity ! to do so; get the number of Channels corresponding ! to diffFreqCentralFreq and multiply by the ! channel width in velocity let veloDiff velo_step*(diffFreqCentralFreq/freq_step) ! add the velocity difference to the central velocity let velo R%HEAD%SPE%VOFF+veloDiff exa velo end procedure freqToVelocity begin procedure show_parts define real part1 part2 let part1 &1 let part2 &2 pen 0 /dashed 3 greg\draw relocate part1 USER_YMIN /user greg\draw line part1 USER_YMAX /user /clip greg\draw relocate part2 USER_YMIN /user greg\draw line part2 USER_YMAX /user /clip pen /def end procedure show_parts !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 ! Construct the line windows from ! ! Input &1 windowSize in MHz !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! begin procedure win ! read in the Line velocities ! for that windows with width 'windowSize' ! are to be created greg1\column x 1 y 2 /file gag_scratch:lines.dat ! lenght of the array storing the window ! number of line (nxy) times 2 for the two window borders ! we have to add +1 since "set wind /var " expects ! the first entry in the array to indicate the number of ! windows def int lenght let lenght nxy*2+1 if (lenght.gt.1) then ! Windows is the array that stores all line Windows if exist(Windows) then del /var Windows endif def real Windows[lenght] /global ! first entry has to be the window number in the array let Windows[1] nxy exa Windows ! Now fill the array with the corresponding window !borders in velocity for i 1 to nxy by 1 say 'x[i]-y[i]' !@freqToVelocity 'x[i]-y[i]' let velo = x[i]-y[i] let Windows['2*i'] velo say 'x[i]+y[i]' !@freqToVelocity 'x[i]+y[i]' let velo = x[i]+y[i] let Windows['2*i+1'] velo next i set wind /var Windows else if (lenght.eq.1) then set wind 0 0 endif end procedure win begin procedure subband set variable spectro read define real jump1 jump2 nfts dum define real ty ty2 ty3 /like ry define integer band reso if exist(velo) then del /var velo endif define int velo /global !define int reso let reso 200 say "Resolution="'reso'"kHz" if ((freq_step.ge.0.194).and.(freq_step.le.(0.200))) then !reso 1 = 200kHz say "Here" exa nchan ! 20727 = number of channels from header exa rchan ! 8961 = reference channel exa fres ! 0.1953 exa vres !if ((fres.gt.0.2).or.(fres.lt.0.19)) then ! pause " ERROR: fres.ne.0.195 MHz " !endif let nfts = (nchan+2)|3 let jump1 = nfts ! edge in number of channels let jump2 = 2*nfts let jump1 = vres*(jump1-rchan)+voff ! edges in velocity let jump2 = vres*(jump2-rchan)+voff else if ((freq_step.ge.0.045).and.(freq_step.le.(0.055))) then exa nchan exa rchan exa fres exa vres if ((fres.gt.0.06).or.(fres.lt.0.04)) then pause " ERROR: fres.ne.50kHz " endif let jump1 = 12494 let jump2 = 12494+12289 ! information from GP, 2-Feb-2012 let jump1 = vres*(jump1-rchan)+voff ! edges in velocity let jump2 = vres*(jump2-rchan)+voff else pause " ERROR: wrong resolution " end if if (jump1.lt.jump2) then continue else let dum = jump1 let jump1 = jump2 let jump2 = dum endif ! this gets the line windows for the current spectrum @win if (debug.eq.1) then set unit v f plot say " Band edges = "'jump1'" "'jump2'" km/s " ! -1227 517 km/s ! define real dum1 dum2 ! let dum1 = jump1 ! let dum2 = jump2 ! set window dum1 dum2 @show_parts 'jump1' 'jump2' @ win draw window say " nchan, rchan, fres, vres "'nchan'" "'rchan'" "'fres'" "'vres' pause "original spectrum with FTS edges" endif ! subtract baselines from each of the three FTS units in the 4GHz band: let ry -1000.00 /where rx.gt.jump1 if (debug.eq.1) then plot pause endif if (debug.ne.1) then base 1 else if (debug.eq.1) then base 1 /plot; draw window; pause "1st FTS unit handled with line window " endif let ty = ry /where rx.lt.jump1 get let ry -1000.00 /where (rx.lt.jump1) let ry -1000.00 /where (rx.gt.jump2) if (debug.eq.1) then plot pause endif if (debug.ne.1) then base 1 else if (debug.eq.1) then base 1 /plot; draw window; pause endif let ty2 = ry /where ((rx.gt.jump1).and.(rx.lt.jump2)) get let ry -1000.00 /where rx.lt.jump2 if (debug.eq.1) then plot pause endif if (debug.ne.1) then base 1 else if (debug.eq.1) then base 1 /plot; draw window; pause "3rd FTS unit handled with line window" endif let ty3 = ry /where rx.gt.jump2 let ry = ty /where rx.lt.jump1 let ry = ty2 /where (rx.gt.jump1).and.(rx.lt.jump2) let ry = ty3 /where rx.gt.jump2 if (debug.eq.1) then plot draw wind !smooth box 10; plot; ! set window dum1 dum2 ! draw window @show_parts 'jump1' 'jump2' @win draw window pause " new spectrum with FTS edges " endif delete /var jump1 jump2 nfts dum ty ty2 ty3 band reso ! delete /var dum1 dum2 end procedure subband !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!