;+
; NAME:
;	wso_write
; PURPOSE:
;	Write magnetic field data array in format of ascii WSO files
; CATEGORY:
;	gen/idl/util
; CALLING SEQUENCE:
	PRO wso_write, data_file, data	, $
		cvariable	= cvariable		, $
		rotation	= rotation		, $
		longitude	= longitude		, $
		header		= header		, $
		fill		= fill			, $
		badfillvalue= badfillvalue	, $
		addfirst	= addfirst		, $
		addlast		= addlast		, $
		hdr_append	= hdr_append	, $
		data_append	= data_append	, $
		obstime		= obstime
; INPUTS:
;	data_file	scalar; type: string
;					
;	data		array[N,30]; type: float
;					N=72 or N=73; magnetic field data
; OPTIONAL INPUT PARAMETERS:
;	rotation=rotation
;					array[N]; type: integer
;						Carrington rotation for each column in Z
;	longitude=longitude
;					array[N]; type: integer
;						heliographic longitude for each column in Z
;	cvariable=cvariable
;					array[N]; type: double
;						Carrington variable = rotation+(1-longitude/360)
;	/addfirst		only used if N=72
;	/addlast		only used if N=72
;						If there are only 72 longitudes we need to add one
;						extra longitude
;						If /addfirst is set then add a copy of the last longitude
;						at the start of the file (with rotation decreased by one).
;						If /addlast is set then add a copy of the first longitude
;						at the end of the file (with rotation increased by one)
;	header=header	scalar or array; type: string
;						is added to the start of the file
;	/fill			before writing to file, bad values (NaNs) are filled using gridfill
;						or is set to badfillvalue (if set).
;						If no filling is done then bad values are written as 'NaN'.
;						(Xuepus program produces isolated NaNs, mostly close to the current sheet.
;						A badfillvalue of zero is the easiest way to get rid of these).
;	badfillvalue=badfillvalue
;					scalar; type: float; default: none
;						used to set bad values (NaNs) is keyword /fill is set
; OUTPUTS:
;	(written to file 'data_file')
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, IsType, gridfill, AngleRange
; SEE ALSO:
;	wso_read
; PROCEDURE:
;	Data for each of the N longitudes are written in four lines using
;	format = '(A2,I4.4,A1,I3.3,8X,6F9.3,3(/,8F9.3))'
;
;	Each group of four records starts with a string CTrrrr:lll
;	where rrrr is the Carrington rotation number (from keyword rotation
;	or integer part of cvariable) and lll the heliographic longitude
;	(from keyword longitude or the fraction of cvariable.
;
;	Longitudes are written from large to small heliographic longitudes
;	(i.e. earlier data first). Latitudes are written from North to South.
;
;	Keyword cvariable takes precedence over keywords rotation and longitude
;	(these keywords will usually be output from wso_read)
; MODIFICATION HISTORY:
;	NOV-2002, Paul Hick (UCSD/CASS)
;	JAN-2003, Paul Hick (UCSD/CASS)
;		Bug fix. The routine replicated the wrong record at the end of the file
;		when the number of longitudes in data is 72.
;		Added /fill keyword.
;	MAY-2003, Paul Hick (UCSD/CASS)
;		Added /addfirst and /addlast keywords. The default is /addlast, so
;		without the keywords the routine behaves as before.
;	DEC-2003, Paul Hick (UCSD/CASS)
;		Added option to write fits file (if extension of output file is
;		.fts, .fit, or .fits). Also allowed appending extra map from keywords
;		hdr_append and data_append. This is used primarily to append the
;		original photospheric synoptic magnetic field map to the source surface
;		map produces by Xuepues hcss program.
;	OCT-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added obstime keyword
;-

InitVar, fill, /key

format = '(A2,I4.4,A1,I3.3,8X,6F9.3,3(/,8F9.3))'

n = size(data, /dim)

IF (n[0] EQ 72 OR n[0] EQ 73) AND n[1] EQ 30 THEN BEGIN

	IF IsType(cvariable, /defined) THEN BEGIN
		rotation  = floor(cvariable)
		longitude = round( (1-(cvariable-rotation))*360 )
		i = where(longitude EQ 0)
		IF i[0] NE -1 THEN BEGIN
			rotation [i] = rotation[i]+1
			longitude[i] = 360
		ENDIF
	ENDIF

	; Check for bad values and fill them if requested.

	IF fill THEN BEGIN
		bad = where(1-finite(data), nbad)
		IF nbad NE 0 THEN BEGIN
			message, /info, strcompress(nbad,/rem)+' bad values filled'
			CASE IsType(badfillvalue, /defined) OF
			0: data = gridfill(data)
			1: data[bad] = badfillvalue
			ENDCASE
		ENDIF
	ENDIF

	extension = strlowcase(GetFileSpec(data_file, part='type'))
	fits_file = ((where(extension EQ ['.fts','.fit','.fits'])) NE -1)[0]

	CASE fits_file OF

	0: BEGIN
		openw, /get_lun, iu, data_file

		FOR i=0,n_elements(header)-1 DO	$			; Add comments at start of file
			printf, iu, header[i]

		; Add magnetic field data in WSO format.
		; Reverse both longitude and latitude order.

		CASE n[0] eq 72 OF
		0: FOR i=n[0]-1,0,-1 DO			$	
				printf, iu, format=format, 'CT',rotation[i],':',longitude[i], reverse(reform(data[i,*]))
		1: BEGIN
			InitVar, addfirst, /key
			InitVar, addlast , /key
			addfirst = addfirst OR 1-addlast

			; If there are only 72 longitudes we need to add one extra longitude
			; If /addfirst is set then add a copy of the last longitude at the start of the file (decrease rotation
			; counter by one). Note that we copy the first column in 'data'. This creates a copy of the last record
			; on the file at the start of the file (with the rotation number decremented by one).
			; If /addlast is set then add a copy of the first longitude at the end of the file (increase rotation
			; counter by one). Note that we copy the last column in 'data'. This creates a copy of the first record
			; on the file at the end of the file (with the rotation number incremented by one).

			IF addfirst THEN BEGIN
				i = 0
				printf, iu, format=format, 'CT',rotation[i]-1,':',longitude[i], reverse(reform(data[i,*]))
			ENDIF

			FOR i=n[0]-1,0,-1 DO			$	
				printf, iu, format=format, 'CT',rotation[i]  ,':',longitude[i], reverse(reform(data[i,*]))

			IF addlast THEN BEGIN
				i = n[0]-1
				printf, iu, format=format, 'CT',rotation[i]+1,':',longitude[i], reverse(reform(data[i,*]))
			ENDIF

		END
		ENDCASE

		free_lun, iu
	END

	1: BEGIN

		data_out = data
		rot_out  = rotation
		lng_out  = longitude

		IF n[0] EQ 72 THEN BEGIN

			InitVar, addfirst, /key
			InitVar, addlast , /key
			addfirst = addfirst or 1-addlast

			IF addfirst THEN BEGIN
				data_out = [data_out, data_out[0,*] ]
				rot_out  = [rot_out , rot_out[0]-1	]
				lng_out  = [lng_out , lng_out[0]	]
			ENDIF

			IF addlast  THEN BEGIN
				i = n[0]-1
				data_out = [data_out[i]  , data_out]
				rot_out  = [rot_out [i]+1, rot_out ]
				lng_out  = [lng_out [i]  , lng_out ]
			ENDIF

		ENDIF

		mkhdr, hdr, IsType(data_out), size(data_out,/dim), extend=IsType(hdr_append,/defined)

		fxaddpar, hdr, 'CARROT'  , rot_out[0], ' Carrington Rotation last column'
		fxaddpar, hdr, 'CARRLONG', lng_out[0], ' Heliographic longitude last column'
		fxaddpar, hdr, 'STEPLONG', round(AngleRange(lng_out[1]-lng_out[0],/pi,/deg)), ' Step size in heliographic longitude'

		i = (where(strpos(header,'OBSERVAT=') EQ 0,complement=j))[0]
		IF i NE -1 THEN BEGIN
			fxaddpar, hdr, 'OBSERVAT', strmid(header[i],strlen('OBSERVAT=')), ' Observatory'
			IF j[0] NE -1 THEN header = header[j] ELSE destroyvar, header
		ENDIF

		i = (where(strpos(header,'OBSTIME=') EQ 0,complement=j))[0]
		IF i NE -1 THEN BEGIN
			obstime = strmid(header[i],strlen('OBSTIME='))
			fxaddpar, hdr, 'OBSTIME', obstime, ' Last time of observation'
			IF j[0] NE -1 THEN header = header[j] ELSE destroyvar, header
		ENDIF

		i = (where(strpos(header,'RSS=') EQ 0,complement=j))[0]
		IF i NE -1 THEN BEGIN
			fxaddpar, hdr, 'RSS', flt_string(header[i]), ' Source surface distance (solar radii)'
			IF j[0] NE -1 THEN header = header[j] ELSE destroyvar, header
		ENDIF

		i = (where(strpos(header,'RCUSP=') EQ 0,complement=j))[0]
		IF i NE -1 THEN BEGIN
			fxaddpar, hdr, 'RCUSP', flt_string(header[i]), ' Cusp distance (solar radii)'
			IF j[0] NE -1 THEN header = header[j] ELSE destroyvar, header
		ENDIF

		i = (where(strpos(header,'BMAP=') EQ 0,complement=j))[0]
		IF i NE -1 THEN BEGIN
			fxaddpar, hdr, 'BMAP', strmid(header[i],strlen('BMAP=')), ' Photospheric map'
			IF j[0] NE -1 THEN header = header[j] ELSE destroyvar, header
		ENDIF

		i = (where(strpos(header,'DESCRIPT=') EQ 0,complement=j))[0]
		IF i NE -1 THEN BEGIN
			fxaddpar, hdr, 'DESCRIPT', strmid(header[i],strlen('DESCRIPT=')), ' Data description'
			IF j[0] NE -1 THEN header = header[j] ELSE destroyvar, header
		ENDIF

		FOR i=0,n_elements(header)-1 DO			$
			IF strtrim(header[i]) NE '' THEN	$
				fxaddpar, hdr, 'COMMENT'+strcompress(i,/rem), header[i]

		writefits, data_file, data_out, hdr

		IF IsType(hdr_append,/defined) THEN $
			writefits, data_file, data_append, hdr_append, /append
	END

	ENDCASE

ENDIF ELSE	$
	message, /info, 'invalid dimensions. No file written.'

RETURN  &  END