pro nsclass,infilelist,obsline,outfilelist,outclassfile ;+ ; NAME: ; NSCLASS ; PURPOSE: ; Converts a list of NEWSTAR generated random-group fits files into generic ; fits files (1 per scan) which can be e.g. imported to CLASS and IRAF. ; Also generates a program file to read in and consolidate all scans in CLASS. ; ; EXPLANATION: ; Fits files written by the SpFITS task in NEWSTAR can be converted into generic ; fits files readable by e.g. IRAF and CLASS. ; ; CALLING SEQUENCE: ; NSCLASS, infilelist, obsline, outfilelist, outclassfile ; ; INPUTS: ; INFILELIST - (string) file name containing the names of the random-group ; fits files produced by NEWSTAR task SpFits. One filename per line. ; OBSLINE - (string) the name of the observed line (transition). ; If not provided, default 'CO(3-2)' is used. ; OUTFILELIST - (string) list of output file names (1 per scan) corresponding ; to all input fits files. If this file does not exist or is empty, ; auto-generated extensions of the input filename are used for the ; output file names. ; OUTCLASSFILE- (string) filename into which to write a CLASS procedure which ; will read in and consolidate the fits files ouput by this IDL ; procedure. If not provided, default name "classread.class" is used. ; ; OUTPUTS: ; 1. generic FITS format files (CLASS/IRAF compatible). 1 file per scan. ; 2. CLASS procedure to read in the fits files output by this IDL ; procedure and consolidate into one master CLASS-format (.bur) ; file for each input NEWSTAR fits file. ; ; ; PROCEDURES CALLED: ; Standard FITS routines in the IDL astronomy library ; ; ; NOTES: ; This is nsclass.pro version 1a, January 2007 ; This routine converts NEWSTAR spectrum fits files to a more portable ; FITS format. The output files have been successfully imported to CLASS and IRAF. ; Use the output CLASS procedure generated by this IDL procedure to read ; in and consolidate the fits files generated by this IDL procedure. ; ; Step-by-step procedure: from e.g. ASTE telescope to CLASS ; 1. The ASTE telescope will produce output data in the raw NEWSTAR format. ; You will probably get a DAT tape with all NEWSTAR fits files. ; e.g. rawdata file name: n1068 ; 2. start NEWSTAR ; 3. in NEWSTAR: read in the DAT tape (use task FROM_CMT with parameter ; "1-99", set the tapedrive name, and execute FROM_CMT). ; 4. in NEWSTAR: use task LINTEG/SPFits to write out each raw NEWSTAR file for ; one filterbank at a time. i.e. specify the input rawdata file from the ; previous step and an output file name, check the box "indiv scan", check ; one filterbank (i.e. A1, A2, A3, A4) at a time, then execute task. ; e.g. you will now have n1068a1.fits n1068a2.fits n1068a3.fits n1068a4.fits ; 5. Use procedure NSCLASS.PRO to read in the NEWSTAR fits files produced ; in the above step, and output CLASS compatible fits files. For help ; on this step, see below. ; 6a. Use the CLASS procedure file generated by NSCLASS.PRO to read in ; and consolidate the fits files produced in the previous step. ; e.g. you will now have one output file, "n1068.bur", which contains ; all scans originally in the raw NEWSTAR file "n1068". ; 6b. start IRAF (change loginuser.cl to include the line ; set imtype =fits ; so that IRAF will work directly with FITS files. ; 7. proceed with processing in CLASS or IRAF. ; ; ; EXAMPLE USAGE OF THIS PROCEDURE: ; IDL> .rnew NSCLASS.PRO ; IDL> NSCLASS,'myinputlist','CO(3-2)','myoutputlist','classread.class' ; ; 'myinputlist' : contains a list of NEWSTAR fits files written with SpFITS. ; one file per line. E.g. contents of file : ; n1068a1.fits ; n1068a2.fits ; 'CO(3-2)' : this is the name of the transition observed. ; 'myoutputlist': is an optional file which contains file names for the output ; fits files. Each input file corresponds to n output files, where n ; is the number of scans in the input fits file. ; E.g. if each input fits file consists of 10 scans, then ; the file 'myoutputlist' should contain 20 lines, as follows: ; n1068a1_0001.fits ; n1068a1_0002.fits ; ... ; n1068a1_0010.fits ; n1068a2_0001.fits ; n1068a2_0002.fits ; ... ; n1068a2_0010.fits ; There is no restriction in the name of the files, as long as you ; have one file per scan. ; The program will then output fits files named as above. ; If file does not exist, or is empty, then default ; extensions are used for the output fits files; these consist of the ; input name followed by "_sc0001" etc. E.g. output files ; n1068a1_sc0001.fits ; n1068a1_sc0002.fits ; ... ; n1068a1_sc0010.fits ; n1068a2_sc0001.fits ; n1068a2_sc0002.fits ; ... ; n1068a2_sc0010.fits ; 'classread.class': this file will contain a class procedure to read in ; and consolidate all fits files produced by NSCLASS.PRO. ; To run it, enter CLASS (part of GILDAS) and at the prompt type: ; LAS> @classread.class ; e.g. you will now have the files n1068a1.bur and n1068a2.bur. ; n1068a1.bur will have all scans originally in the raw NEWSTAR ; file n1068a1, and the same applies to n1068a2. ; ; REVISION HISTORY: ; Juan Cortes, Universidad de Chile. Created in April 2006 ; Laura Perez, fixed bugs. April 2006 ; Documentation update J. Cortes Jun 2006 ;- ; N. Nagar, UdeC, revised parameters and program control, Dec 2006 ;- ; procedure starts here --------- print,'nsclass,infilelist,obsline,outfilelist,outclassfile' ; check input file exists; exit if it doesn't s1=0 & nlines='' openr,1,infilelist while (not eof(1)) do begin readf,1,nlines s1=s1+1 ; counts number of input NEWSTAR fits files endwhile close,1 if (s1 EQ 0) then begin print, 'Empty . Please try again" return endif infile=strarr(s1) ; names of all input NEWSTAR fits files openr,1,infilelist & readf,1,infile & close,1 ; check if output fits file list exists; if not, use default names s2=0 ; count for number of generic output FITS files res=findfile(outfilelist,count=fcount) if (fcount GT 0) then begin ; found a file, check if nonempty nlines='' openr,1,outfilelist while (not eof(1)) do begin readf,1,nlines if (strtrim(nlines,2) NE '') then s2=s2+1 ; increment if line is nonblank endwhile close,1 endif ; if fcount GT 0 if ((fcount EQ 0) OR (s2 EQ 0)) then begin ; no file, or empty/blank file print, 'Could not find output fits name list - Will use default file naming' endif else begin ; use the output file naming from file outfile=strarr(s2) ; names of all output files openr,1,outfilelist & readf,1,outfile & close,1 endelse ; check existence of transition, else use default if (obsline EQ '') then obsline='CO(3-2)' ; check existence of name, and then open, the class program file if (outclassfile EQ '') then begin print,"No classfile provided, using default name classread.class" outclassfile = 'classread.class' endif openw,5,outclassfile vlight=2.99792458e+08 ; vel of light in meter/sec totcount=0 ; initialize running count for number of fits files ouput ; header comments in outclassfile printf,5,'! to run this file, enter CLASS and type the following line' printf,5,'! @',outclassfile ; read in fits files, convert header formats, and write out within this loop for filecnt=0,(size(infile))(1)-1 do begin ;for each input NEWSTAR fits file in the filelist ; define name of final CLASS .bur file outclassname=infile(filecnt) pos = strpos(outclassname,".fits") pos2= strpos(outclassname,".FITS") if (pos GE 0) then $ outclassname= strmid(outclassname,0,pos) $ else begin if (pos2 GE 0) then outclassname= strmid(outclassname,0,pos2) endelse printf,5,'file out ',outclassname," new" print,'----- processing file: ',infile(filecnt),'-----' ; Notes on NEWSTAR fits file headers, useful for future updates of this .pro: ;1. in NEWSTAR fits, a few header parameter names and value pairs are listed ; explicitly together, e.g. CRPIX, NAXIS2, BSCALE. ; Other header parameters names are listed in the header without their values. ; For these, the corresponding values are stored in each scan's data array, since ; the values are often different for each scan within the NEWSTAR fits file. ; Here, for each input NEWSTAR fits file, the data array is read into a variable ; named "data". "data" is a two dimensional array: the first dimension corresponds ; to the scan number, e.g. data[1] corresponds to the 1st scan. ; data[i].params[j] will give the j'th header parameter value of the i'th scan. ; data[i].array[j] will give the j'th data value of the i'th scan. ;2. note that the NEWSTAR fits file has freq/vel as axis_2. We need to output ; freq/vel as axis_1. ;3. The NEWSTAR fits header contains a rest frequency and a tracking frequency. ; the tracking frequency is the same for all arrays A1 to A4. It is the ; frequency of the receiver LO (?). ; The rest frequency refers to the center frequency of each array A1 to A4. ; if the arrays were offset w.r.t. each other, then the rest frequencies of ; each array will be different. ;4. In NEWSTAR fits file: if CTYPE2=FREQ then CRVAL2 in the NEWSTAR ; fits will be the same as the rest frequency (Hz). CDELT2 in the NEWSTAR ; header will be freq per channel in Hz (at ASTE this is 500kHz/ or 125kHz/channel ; depending on whether the highBW or lowBW mode was used in the backend). ; In NEWSTAR fits file: if CTYPE2=VELO-LSR then CRVAL2 in the NEWSTAR fits will ; be the VLSR (or Vhelio) used in meter/sec. CDELT2 will be the velocity ; per channel in meter/sec. ;5. data[0].array[1023], i.e. last data value of the 1st scan in a NEWSTAR fits ; file is always = -2147483648 (NEWSTAR value for "BLANK"). While the value ; of "BLANK" is specified in the ouput header, several softwares do not ; recognize and correct this header parameter. So, for now, set this value ; to that of data[0].array[1022]. Look for a better fix for this later. ; ; Notes on CLASS fits headers, useful for future updates of this .pro: ; in CLASS, the frequency of a channel is given by ; F(i) = RESTFREQ + CRVAL1 + ( i - CRPIX1 ) * CDELT1 ; so we need to set CRVAL1 = 0, and put our restfreq value in RESTFREQ ; velocity of a channel: V(i) = VLSR + ( i - CRPIX1 ) * DELTAV ; so we need to set VLSR, and DELTAV even if the channel axis is frequency ; Notes on IRAF fits headers: ; the frequency of each channel is given by ; F(i) = CRVAL1 + ( i - CRPIX1 ) * CDELT1 ; i.e. different from class. Thus we either have to live with a strange ; channel axis in IRAF, or add an input to nsclass to specify whether the ; user wants CLASS- or IRAF-compatible output headers. data=mrdfits(infile(filecnt),0,h1) ; read in fits file help,/struct,data ; print details of "data" to screen gcount=sxpar(h1,'GCOUNT') ; number of scans within the FITS file pcount=sxpar(h1,'PCOUNT') ; number of header parameter values for each scan pscale=dblarr(pcount) ; scaling factor for each header parameter value crpix=sxpar(h1,'CRPIX2') nel=sxpar(h1,'NAXIS2') ; number of frequency/velocity channels. bscale=sxpar(h1,'BSCALE') ; generate the explicit NEWSTAR header names used for the scaling parameters. ; and read in the scaling parameter values from the header for i=1,pcount do begin string_i=strtrim(i,2) pscale(i-1)=sxpar(h1,'PSCAL'+string_i) endfor value=dblarr(gcount,pcount) ; array to hold the header parameter values ; read in header parameter values and scale them by their resp. scaling factor for i=1,pcount do begin for j=0,gcount-1 do begin value(j,i-1)=data[j].params[i-1]*pscale[i-1] endfor endfor cube=lonarr(nel,gcount,1) ; array to hold the data values of all scans help,cube history=sxpar(h1,'history') nelhist=n_elements(history) for i=0,gcount-1 do begin ; for each scan in a given input fits file for j=0,nel-1 do begin cube(j,i,0)=long(data[i].array[j]) ; read this scan's data values endfor ; temp fix for the problem that data[0].array[1023] always = BLANK if ((i EQ 0) AND (cube(1023,i,0) EQ -2147483648)) then $ cube(1023,i,0)=cube(1022,i,0) ; header values for the current scan of the NEWSTAR fits file crval2=value[i,4] cdelt2=value[i,6] crpix2=0. crval3=value[i,5] cdelt3=value[i,7] crpix3=0. bzero=sxpar(h1,'BZERO') bunit=sxpar(h1,'BUNIT') blank=sxpar(h1,'BLANK') ctype2=sxpar(h1,'CTYPE2') object=sxpar(h1,'OBJECT') equinox=sxpar(h1,'EQUINOX') raoj=value[i,4] decoj=value[i,5] dra=value[i,6] ddec=value[i,7] cubetemp=lonarr(nel,10,10) fxhmake,headertemp,cubetemp header=headertemp ; this is the header to be written out date=value[i,15] ; fix data format in the following commands idate=long(date) year=long(date/10000) fdate=idate-year*10000 month=long(fdate/100) day=fdate-month*100 year1=strtrim(year,2) month1=strtrim(month,2) day1=strtrim(day,2) if (year lt 10) then year1='0'+year1 if (year lt 80) then year1='20'+year1 else year1='19'+year1 if (month lt 10) then month1='0'+month1 if (day lt 10) then day1='0'+day finaldate=year1+'-'+month1+'-'+day1 ; set values in the header to be written out sxaddpar,header,'naxis',3 sxaddpar,header,'naxis2',1 sxaddpar,header,'naxis3',1 sxaddpar,header,'blank',blank sxaddpar,header,'bscale',bscale sxaddpar,header,'bzero',bzero sxaddpar,header,'bunit',bunit sxaddpar,header,'scan-num',i+1 fxaddpar,header,'crpix1',crpix ; earlier version (JC?) had the foll definition for CDELT1 delfreq=-1.0*(value[i,1]/vlight)*(value[i,12]/(1.0+(value[i,1]/vlight))) ; this is probably -1 * (del_v/c) * (v / (1 + del_v/c)) ; or with nu instead of v, which would make it more like the optical f_sky ; we don't use this value anymore ; some complexity in setting CRVAL1 and CRDELT1 i.e. freq or vel axis if (strtrim(ctype2,2) EQ 'FREQ') then begin ; channel axis is in FREQuency fxaddpar,header,'ctype1',ctype2 ;fxaddpar,header,'crval1',value[i,0],format='g16.9' ; need this for IRAF fxaddpar,header,'crval1',0.,format='g16.9' ; crval=0 for CLASS, see notes above fxaddpar,header,'cdelt1',value[i,1],format='g16.9' ; If crtype1=freq: we need to calculate the delta velocity of each channel, ; DELTAV, to keep CLASS happy. ; Since the ASTE bandwidth (128 or 512MHz per array) is much smaller ; than nu ( ~ 350GHz), we can use del_nu/nu = del_v/c delvel= -1.0 * vlight * value[i,1]/value[i,12] ; value in meter/sec ; have to check if setting VLSR = 0 is always valid - NN. fxaddpar,header,'VLSR',0. fxaddpar,header,'DELTAV',delvel,'VELOCITY PER CHANNEL: VALUE IS METER/SEC',format='g16.9' endif ; NN needs to re-check this part with more datasets (in VELO-LSR) ; e.g. if VELO-LSR is used, are arrays shifted w.r.t. each other by using different ; values of VELO-LSR. Or is the shift reflected by different RESTFREQs? ; please send updates/comments to NN ; for now, the axis is always ouput as CTYPE=FREQ if (strtrim(ctype2,2) NE 'FREQ') then begin ; channel axis is in VELO-LSR (or VELO-HEL?) fxaddpar,header,'ctype1','FREQ' ; change channel axis from VELO to FREQ ; f_sky is from optical defn: = f_rest * (1 / (1+v/c)) f_sky= value[i,12] * (1.0 / (1.0 + value[i,0]/vlight)) f_offset= f_sky - value[i,12] ; set crval1 to freq offset implied by VLSR (works for CLASS) fxaddpar,header,'crval1',f_offset,'FREQ OFFSET (Hz) TO REFLECT VLSR USED',format='g16.9' ; in IRAF, crval1 should be f_sky, as follows ; fxaddpar,header,'crval1',f_sky,'SKY FREQ',format='g16.9' ; del_freq/freq = del_v / c delfreq= -1.0 * f_sky * (value[i,1]/vlight) fxaddpar,header,'cdelt1',delfreq,'FREQ PER CHANNEL: VALUE IN HZ',format='g16.9' fxaddpar,header,'VLSR',value[i,0],'V_LSR: VALUE IN METER/SEC',format='g16.9' fxaddpar,header,'DELTAV',value[i,1],'VELOCITY PER CHANNEL: VALUE IN METER/SEC',format='g16.9' endif fxaddpar,header,'ctype2','RA' ;'RA---ARC';'RA--TAN' fxaddpar,header,'crpix2',crpix2,format='g16.9' fxaddpar,header,'crval2',crval2,format='g16.9' fxaddpar,header,'cdelt2',cdelt2,format='g16.9' fxaddpar,header,'ctype3','DEC'; 'DEC--ARC';'DEC---TAN' fxaddpar,header,'crpix3',crpix3,format='g16.9' fxaddpar,header,'crval3',crval3,format='g16.9' fxaddpar,header,'cdelt3',cdelt3,format='g16.9' fxaddpar,header,'object',object fxaddpar,header,'equinox',equinox,format='g16.9' fxaddpar,header,'RA',raoj,'RA OF MAP CENTER: VALUE IS DEGREE',format='g16.9' fxaddpar,header,'DEC',decoj,'DEC OF MAP CENTER: VALUE IS DEGREE',format='g16.9' fxaddpar,header,'RAOFF',dra,'RA OFFSET: VALUE IS DEGREE',format='g16.9' fxaddpar,header,'DECOFF',ddec,'DEC OFFSET: VALUE IS DEGREE',format='g16.9' fxaddpar,header,'obstime',value[i,17],'INTEGRATION TIME: VALUE IS SEC',format='g16.9' fxaddpar,header,'navsp',value[i,18],'NO OF SPECTRA AVERAGED',format='g16.9' fxaddpar,header,'ARRAY-NO',value[i,16],'AOS ARRAY NUMBER',format='g16.9' fxaddpar,header,'RESTFREQ',value[i,12], 'CENTER FREQUENCY: VALUE IS HZ',format='g16.9' fxaddpar,header,'FQTRK',value[i,13],'TRACKING FREQUENCY: VALUE IS HZ',format='g16.9' fxaddpar,header,'TAU',value[i,24],'SKY OPTICAL DEPTH' fxaddpar,header,'DATE-OBS',finaldate fxaddpar,header,'LINE',obsline ; These values are in the NEWSTAR FITS header. Pass them on to the FITS header. ; these are the absolute RA and DEC of the current offset position fxaddpar,header,'RAOFFSET',value[i,8],'RA OF OFFSET POSITION: VALUE IS DEGREE',format='g16.9' fxaddpar,header,'DECOFFSET',value[i,9],'DEC OF OFFSET POSITION: VALUE IS DEGREE',format='g16.9' ; add history to output header for histcnt=0,nelhist-1 do begin fxaddpar,header,'HISTORY',history(histcnt) endfor fxaddpar,header,'HISTORY','Converted to generic FITS by NSCLASS.PRO (Ver 1a, Jan2006)' if (strtrim(ctype2,2) NE 'FREQ') then $ ; channel axis was not in FREQ fxaddpar,header,'HISTORY','NSCLASS.PRO: Original CTYPE=VELO-LSR; has been changed to FREQ' ; output fits file: naming and writing if ((s2 NE 0) AND (totcount LT s2)) then begin ; output names provided outfilename = outfile(totcount) endif else begin ; no output names provided outfilename=infile(filecnt) pos = strpos(outfilename,".fits") pos2= strpos(outfilename,".FITS") if (pos GE 0) then $ outfilename= strmid(outfilename,0,pos) $ else begin if (pos2 GE 0) then outfilename= strmid(outfilename,0,pos2) endelse outfilename= outfilename + '_sc' + string((i+1),format='(I4.4)') + '.fits' endelse if ((s2 NE 0) AND (totcount GE s2)) then $ print,"run out of output file names. will use default= ",outfilename writefits,outfilename,long(cube(*,i,0)),header printf,5," las\fits read ",outfilename printf,5," write" totcount=totcount+1 endfor ; for i (each scan of a given input fits file) printf,5," ! end of file", infile(filecnt) stop endfor ; for filecnt (each input fits file) printf,5,"exit" close,5 print,"-----------------------------------------------------------" print,"NEWSTAR to generic FITS conversion completed successfully: " print,"Number of input fits files: ",filecnt print,"Number of output fits files: ",totcount print,"Use ",outclassfile," to read and consolidate output fits files in CLASS" print," " print,"Output channel axis is CRTYPE1=FREQ " print,"If you display these spectra with CLASS: ---" print," If your original NEWSTAR file had channel axis CRTYPE=FREQ " print," then reference channel (CRPIX) is set to VEL=0 and FREQ=your sky freq" print," If your original NEWSTAR file had channel axis CRTYPE=VELO-LSR:" print," then reference channel (CRPIX) is set to Velocity= V_LSR and FREQ=F_sky" print,"If you display these spectra with IRAF: ---" print," If your original NEWSTAR file had channel axis CRTYPE=FREQ" print," then reference channel (CRPIX) is set to Frequency=0" print," If your original NEWSTAR file had channel axis CRTYPE=VELO-LSR" print," then reference channel (CRPIX) is set to Frequency= F_offset" print," " print,"here: F_offset = F_sky - RESTFREQ and F_sky = RESTFREQ * (1/(1+VLSR/c))" print," " print,"Suggestions/bugs: jcortes@das.uchile.cl, nagar@astro-udec.cl" print,"-----------------------------------------------------------" end