;+ ; NAME: ; Radio Astronomy Group Fits Write ; PURPOSE: ; Writes RAG spectrograms in FITS files ; including non-regular axes. ; CALLING SEQUENCE: ; RAGFitsWrite, image [, xAxis, yAxis ] ; INPUTS: ; image: 2D or 3D array of any type except string and complex. ; xAxis, yAxis: 1D arrays of any type but string and complex. ; They must have the same dimensions as their ; corresponding dimension in "image". ; The third dimension of the image must have 2 elements: ; (0) the flux and (1) the polarization ; KEYWORDS: ; FILENAME: string containing name of the file with ; extension ; CONTENT: default: 'SPECTROGRAM' ; ORIGIN: default; 'RAG / ETH Zurich' ; TELESCOPE: default: 'Bleien Radio' ; INSTRUMENT: default: 'PHOENIX' ; OBJECT: default: 'SUN' ; DATEOBS:observation date in form dd/mm/yy, ; default: empty string. ; TIMEOBS:start time of observation ( hh:mm:ss.ddd ) ; default: the first element of the x-axis ; This keyword should generally not be used, ; since the time start is generally contained in ; first element of the x-axis. If timeobs is used, ; its quantity will be added to the x(0)-element ; when the file is read by RagFitsRead. ; DATEEND, TIMEEND: as DATEOBS and TIMEOBS, ; but for the end of the observation. ; default: TIMEEND = TIMEOBS; ; TIMEEND = empty string. ; BZERO, BSCALE: to set if the data has to be reconstructed ; with data = BZERO + BSCALE*image ; Default: BZERO = 0; BSCALE = 1 (see ; the keyword SCALING. ; XZERO, YZERO: as BZERO but for x, resp, y axes. ; YSCALE, YSCALE: as BSCALE but for x, resp. y axes. ; These 4 keywords are only used in case of non- ; regular axes. ; SCALING: If present, the procedure scales automatically ; the data and writes in the fits header ; the computed BSCALE AND BZERO values. ; If not present, the current type of the array is ; stored. ; "Byte": the data is written in 8-bit format. ; "Integer": the data is written in 16-bit format. ; "LongInt": the data is written in 32-bit format. ; BUNIT: string containing the units of the image pixels. ; default: '45*LOG(SFU+10)' ; CTYPE1, CTYPE2: string containing the FITS type of ; physical coordinates. Defaults: ; CTYPE1 = 'TIME' and CTYPE2 = 'FREQ' ; COMMENT, HISTORY: strings for description ; of the image ; INTERACTIVE: if set, the filename is asked interactively ; BYTE: If SCALED is not used, BYTE allows to write ; real data in byte format (symply transforming ; the data without byte scaling ) ; NOROUND : if set, the axis is NOT rounded to two significant ; digits, (default for RAG fits ) ; SIDE EFFECT: ; A fits file is created. ; MODIFICATION HISTORY ; Created: A. Csillaghy, ETHZ, October 1992 ; BYTE keyword added in May 93, A.Cs. ; NOROUND in Nov 93, A.Cs ;- PRO AxisKeywords, kwrdNb, axis, nEls, typeName, $ primaryHeader, regular, cDelt, NOROUND = noRound br = FindBreaks( axis,/NOTREG ) IF br(0) EQ -1 THEN regular = 0B ELSE regular = 1B crVal = axis(0) crPix = 0 IF nEls GT 1 THEN BEGIN IF regular NE 1 THEN BEGIN IF kwrdNb EQ 1 THEN axisName = 'time' $ ELSE axisName = 'frequency' FxAddPar, primaryHeader, 'COMMENT', $ 'The '+ axisName + ' axis is not regular' ENDIF cDelt = axis(1) - axis(0) IF cDelt NE 0 AND NOT Keyword_Set( NOROUND ) THEN BEGIN i = 0 ; Round the step to 2 significant numbers WHILE (Abs(Fix( cDelt )) LT 10) DO BEGIN i = i+1 cDelt = cDelt*10 ENDWHILE IF i GT 0 THEN BEGIN cDelt = Round( cDelt ) / 10.0^i FxAddPar, primaryHeader, 'COMMENT', $ 'WARNING: the value of CDELT' + kwrdNb + $ ' may be rounded ! ' ENDIF ENDIF ENDIF ELSE cDelt = 0 FxAddPar, primaryHeader, "CRVAL"+ kwrdNb, crVal, $ ' value on axis ' + kwrdNb + ' at the reference pixel' FxAddPar, primaryHeader, "CRPIX"+ kwrdNb, crPix, $ ' reference pixel of axis ' + kwrdNb FxAddPar, primaryHeader, "CTYPE"+ kwrdNb, typeName, $ ' title of axis ' + kwrdNb FxAddPar, primaryHeader, "CDELT"+ kwrdNb, cDelt, $ ' step between first and second elements in axis' END ; write axis information ;---------------------------------------------- FUNCTION PackAxis, zero, scale, axis IF (N_Elements( zero ) NE 0) AND (N_Elements( scale ) NE 0) $ THEN RETURN, Long( Double( axis - zero ) / scale )$ ELSE RETURN, Pack( axis, zero, scale, /LONGINT ) END; pack axis ; ---------------------------------------------- PRO MakeBinTable, xAxis, yAxis, nx, ny, $ fileName, xZero, yZero, xScale, yScale IF N_Elements( xZero ) NE 0 THEN xZero = Dbl(xZero) IF N_Elements( yZero ) NE 0 THEN yZero = Dbl(yZero) IF N_Elements( xScale ) NE 0 THEN xScale = Dbl(xScale) IF N_Elements( yScale ) NE 0 THEN yScale = Dbl(yScale) FxBHMake, header, 1, 'Axes', /INIT xAxisPacked = PackAxis( xZero, xScale, xAxis ) yAxisPacked = PackAxis( yZero, yScale, yAxis ) FxBAddCol, indexX, header, xAxisPacked, 'time', TUNIT = 's', $ TDMIN = Min( xAxis ) , TDMAX = Max( xAxis ), $ TZERO = xZero, TSCAL = xScale FxBAddCol, indexY, header, yAxisPacked, 'frequency', $ TUNIT = 'MHz', $ TDMIN = Min( yAxis ) , TDMAX = Max( yAxis ), $ TZERO = yZero, TSCAL = yScale FxBCreate, unit, fileName, header FxBWrite, unit, xAxisPacked, indexX, 1 FxBWrite, unit, yAxisPacked, indexY, 1 FxBFinish, unit END; make binary table ; ----------------------------------------------- PRO RAGFitsWrite, image, xAxisP, yAxisP, FILENAME = filename, $ CONTENT = content, ORIGIN = origin,$ TELESCOPE = telescope, INSTRUMENT = instrument, $ OBJECT= object, BYTE = byte, $ DATEOBS = dateObs, TIMEOBS = timeObs, $ DATEEND = dateEnd, TIMEEND = timeEnd, $ BZERO = bZero, BSCALE = bScale, BUNIT = bUnit, $ CTYPE1 = cType1, CTYPE2 = cType2, $ COMMENT = comment, SCALING = scaling, $ HISTORY = history, INTERACTIVE = interactive, $ XZERO = zXero, YZERO = yZero, $ XSCALE = xScale, YSCALE = yScale, NOROUND = noRound ; On_Error, 2 ; 1. Parameter and keyword tests IF NOT Keyword_Set( FILENAME ) OR $ Keyword_Set( INTERACTIVE )THEN BEGIN Print, 'FITS File Write' filename = '' Read, 'Please enter the name of the file: ', fileName ENDIF fnSize = Size( fileName ) IF fnSize(0) NE 0 OR fnSize(1) NE 7 THEN $ Message, 'File name must be a scalar string ', /INFO IF N_Elements( image ) EQ 0 THEN $ Message, 'An image must be defined', /INFO nx = N_Elements( image(*,0,0)) ny = N_Elements( image(0,*,0)) nz = N_Elements( image(0,0,*)) IF N_Elements( xAxisP ) EQ 0 THEN xAxis = IndGen( nx ) $ ELSE xAxis = xAxisP IF N_Elements( yAxisP ) EQ 0 THEN yAxis = IndGen( ny ) $ ELSE yAxis = yAxisP minImage = Min( image, MAX = maxImage ) IF NOT Keyword_Set( CONTENT ) THEN $ content = 'Spectrogram' IF NOT Keyword_Set( ORIGIN ) THEN $ origin = 'RAG / ETH Zurich' IF NOT Keyword_Set( TELESCOPE ) THEN $ telescope = 'Bleien Radio' IF NOT Keyword_Set( INSTRUMENT ) THEN $ instrument = 'Phoenix' IF NOT Keyword_Set( OBJECT ) THEN object = 'Sun' IF NOT Keyword_Set( DATEOBS ) THEN dateObs = '' maxTime = HMSConvert( Max(xAxis) ) IF NOT Keyword_Set( TIMEOBS ) THEN BEGIN timeObs = HMSConvert( xAxis(0) ) xAxis = xAxis - xAxis(0) ENDIF IF NOT Keyword_Set( DATEEND ) THEN dateEnd = dateObs IF NOT Keyword_Set( TIMEEND ) THEN BEGIN timeEnd = maxTime ENDIF IF NOT Keyword_Set( BZERO ) THEN bZero = 0 IF NOT Keyword_Set( BSCALE ) THEN bScale = 1 IF NOT Keyword_Set( BUNIT ) THEN bUnit = '45*LOG(SFU+10)' IF NOT Keyword_Set( CTYPE1 ) THEN cType1 = 'TIME' IF NOT Keyword_Set( CTYPE2 ) THEN cType2 = 'FREQ' IF NOT Keyword_Set( CTYPE3 ) THEN cType3 = 'MODE' IF NOT Keyword_Set( COMMENT ) THEN comment = '' IF NOT Keyword_Set( HISTORY ) THEN history = '' ; 2. write primary header and array Print, 'Writing FITS file ... ' ; scale the array scaled = Keyword_Set( SCALING ) IF scaled THEN BEGIN type = Size( scaling ) IF type(1) NE 7 THEN scaling = 'I' bigData = 0; N_Elements( image ) GT 500000L with new memory not necessary CASE StrUpcase( StrMid( scaling, 0, 1 ) ) OF "B": image = Pack( image, bZero, bScale, /BYTE, $ BIGDATA = bigData ) "L": image = Pack( image, bZero, bScale, /LONGINT, $ BIGDATA = bigData ) ELSE: image = Pack( image, bZero, bScale, /INTEGER, $ BIGDATA = bigData ) ENDCASE ENDIF ELSE IF Keyword_Set( BYTE ) THEN $ image = Byte( image ) ; create primary header FxHMake, primaryHeader, image, /EXTEND, /DATE, /INIT FxAddPar, primaryHeader, 'CONTENT', content, $ 'Title of image' FxAddPar, primaryHeader, 'ORIGIN', origin, $ 'Organization name' FxAddPar, primaryHeader, 'TELESCOP', telescope, $ 'Location of the telescope' FxAddPar, primaryHeader, 'INSTRUME', instrument, $ 'Name of the spectrometer' FxAddPar, primaryHeader, 'OBJECT', object FxAddPar, primaryHeader, 'DATE-OBS', dateObs, $ ' date observation starts' FxAddPar, primaryHeader, 'TIME-OBS', timeObs, $ ' time observation starts' IF Keyword_Set( DATEEND ) THEN $ FxAddPar, primaryHeader, 'DATE-END', dateEnd, $ ' date observation ends' FxAddPar, primaryHeader, 'TIME-END', timeEnd, $ ' time observation ends' FxAddPar, primaryHeader, 'BZERO', bZero, ' scaling offset' FxAddPar, primaryHeader, 'BSCALE', bScale, 'scaling factor' FxAddPar, primaryHeader, 'BUNIT', bUnit, 'z-axis title' FxAddPar, primaryHeader, 'DATAMIN', Double(minImage ), $ 'Minimum element in image' FxAddPar, primaryHeader, 'DATAMAX', Double(maxImage), $ 'Maximum element in image' AxisKeywords, '1', xAxis, nx, cType1, primaryHeader, xRegular, dx, $ NOROUND = noRound AxisKeywords, '2', yAxis, ny, cType2, primaryHeader, yRegular, dy, $ NOROUND = noRound nLines = N_Elements( comment ) FOR i=0,nLines-1 DO $ FxAddPar, primaryHeader, 'COMMENT', comment(i) nLines = N_Elements( history ) FOR i=0,nLines-1 DO $ FxAddPar, primaryHeader, 'HISTORY', history(i) IF scaled THEN $ FxWrite, fileName, primaryHeader, image, /NOUPDATE $ ELSE FxWrite, fileName, primaryHeader, image ; 2 write binary table for axes IF (xRegular NE 1) OR (yRegular NE 1) THEN $ MakeBinTable, xAxis, yAxis, nx, ny, $ fileName, xZero, yZero, xScale, yScale END ; fits write