pro YOHKOH, data=data, bin=RREBIN,justbin=JustBin, nocalc=NoCalc

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	YOHKOH
; PURPOSE:
;	Convert the IDL save file SYNMAPS.SAV to ASCII files
; CATEGORY:
;	I/O
; CALLING SEQUENCE:
;	yohkoh, bin=4
; INPUTS:
;	IDL save file SYNMAPS.SAV
;	If its not in $LIB, the file is on the backup optical disk HOSA74-1
; OPTIONAL INPUT PARAMETERS:
;	REBIN		rebin factor (valid values are 1,2,4,8,16; default is 4)
; OUTPUTS:
;	ASCII files in working directory
;	The output files names are ????.YOH, where ???? is the Carrington
;	rotation number
; CALLS:
;	CarringtonT, NewcombSun, BZero, gridgen
; PROCEDURE:
; >	The SYNMAPS.SAV file is restored using the IDL RESTORE command.
;	The array DATA(624,256,16) is written into 16 separate ASCII files,
;	one file for each synoptic map
; >	The output files can be read and plotted by the Fortran SANYOH.EXE
;	program. The current SANYOH version assumes the default rebin value of
;	4; if another value is used the program must be adapted accordingly,
;	i.e. change the parameter declarations for MAXLNG and MAXLAT.
; >	The original maps have dimensions 624 by 256. The default REBIN value
;	reduces them to 156 by 64 before writing them to ASCII files.
; MODIFICATION HISTORY:
;	JUN-1993, Paul Hick (UCSD)
;-

ArcsecPerPix = 9.812			; # arcseconds per pixel

if not keyword_set(RREBIN) then RREBIN = 16
message, /info, 'Rebinning value is '+strcompress(RREBIN,/rem)

A = where ( [1,2,4,8,16] eq RREBIN )
if A[0] eq -1 then message, 'Invalid rebin value. Valid values are 1,2,4,8,16

sdata = size(data)

Nlng  = sdata[1]
Nlat  = sdata[2]
Nmaps = sdata[3]

help, data
print, "  # dimensions : ",sdata[0]
print, "  # maps       : ",Nmaps

if keyword_set(NoCalc) then return

CRstart = 1847   &  echo, 'Start Carrington rotation', CRstart
CRend = CRstart+Nmaps-1
CRNUM = CRstart+indgen(Nmaps)

if n_elements(CRstart) eq 0 then CRstart = CRNUM[0]
CRstart = ((CRstart-CRNUM[0]) > 0) < (Nmaps-1)

if n_elements(CRend) eq 0 then CRend = CRNUM[Nmaps-1]
CRend = ((CRend-CRNUM[0]) > 0) < (Nmaps-1)

if not keyword_set(JustBin) then begin
	xlng = gridgen(Nlng, range=[0,360])		; Heliographic long
	ylat = gridgen(Nlat, range=90*[-1,1])	; Heliographic lat
endif

for i=CRstart,CRend do begin		; Loop over Carrington rotations

	if not keyword_set(JustBin) then begin

; Get times TT at which Earth was at longitude xlng in rotation CRNUM(i)
; Then use TT to calculate B0 angle and solar radius

		TT = CarringtonT(CRNUM[i]+(1-xlng/360))
		B0 = BZero(TT, /degrees)			; B0 angle (in degrees)
		j  = NewcombSun(TT, /distance, diskradius=R0)
		R0 = R0/ArcsecPerPix					; Radius of Sun (in pixels)

		message, /info, 'Warping Carrington rotation '+strcompress(CRNUM[i],/rem)+' ...'
		message, /info, 'Average Bzero = '+string(total(B0)/n_elements(B0))

; Linear interpolation to go from equal steps in sin(lat) to equal steps in latitude

		for j=0,Nlng-1 do begin			; Loop over all longitudes

; Create a grid that is equally spaced in sin(Heliographic lat)

			ysin = ( findgen(Nlat)-(Nlat-1)/2. ) / R0[j]

; Assume that Greg shifted each strip over a distance sin(B0) (south for
; neg B0, north for pos B0). Take this shift into account by adjusting the
; 'yold' array: the center of the disk, i.e. the center of the yold array
; has value -sin(B0) (the center of ysin=0)

			yold = ysin-sin(B0[j]/!radeg)
			vold = data[j,*,i]			; Extract fnc-values

; Drop the part where abs(yold) gt 1. (The intensities in the part that is
; dropped is zero anyway).

			valid = where(abs(yold) le 1.)
			yold = yold[valid]  &  vold = vold[valid]
			yoldmin = min(yold)  &  yoldmax = max(yold)

; Find the fnc-values in a grid that is equally space in latitude (degrees)

			ynew = ylat-B0[j]
			snew = sin(ynew/!radeg)

			if B0[j] ge 0 then begin
				nvalid = where (snew gt yoldmax)
				if invalid[0] ne -1 then		$
					message, /info, strcompress(CRNUM[i],/rem)+'-'+	$
						strcompress(xlng[j],/rem)+' , B0='+		$
						strcompress(B0[j],/rem)+' : lost '+		$
						strcompress(B0[j]+asin(yoldmax)*!radeg,/rem)+	$
						' to 90 deg'

			endif else begin
				invalid = where (snew lt yoldmin)
				if invalid[0] ne -1 then		$
					message, /info, strcompress(CRNUM[i],/rem)+'-'+	$
						strcompress(xlng[j],/rem)+' , B0='+		$
						strcompress(B0[j],/rem)+' : lost -90 to '+	$
						strcompress(B0[j]+asin(yoldmin)*!radeg,/rem)+	$
						' deg'

			endelse

			vnew = interpol(vold, yold, snew)
			data[j,*,i] = round(vnew)

		endfor

	endif else	$
		message, /info, 'Just binning rotation '+strcompress(CRNUM[i],/rem)

	map = data[*,*,i]
	if RREBIN ne 1 then map = rebin ( map , Nlng/RREBIN, Nlat/RREBIN)
	mdata = size( map )

; Open fixed record length ASCII file; write each map to separate file.
; Each longitude on a seperate record : mdata[1] records
; Each record contains mdata[2] latitudes

; File name is Carrington rotation number

	RECL = 4*mdata[2]
	FORMT = '('+strcompress(mdata[2],/remove_all)+'I4)'

	file = string(CRNUM[i],format='(i4.4)')+'.YOH;'
	message, /info, 'Creating file  '+file
	openw, unit, file, /get_lun, RECL, width=RECL
	for j=mdata[1]-1,0,-1 do printf, unit, format=FORMT, map[j,*]
	printf, unit, '_______'			; Write array and add last record
	free_lun, unit				; Close file

endfor

return  &  end