;+
; ZHELIOGRAPHINIZE
;
; PURPOSE - To remap an image from image to heliographic co-ordinates
;
; USAGE
;	zheliographinize,in=in,ihdr=ihdr,out=out,ohdr=ohdr,soho=soho $
;		,date=date,t0=t0,[/howard],[/rigid]
;
; REQUIRED KEYWORDS:
;	IN,IHDR - the input image and header to heliographinize
;	OUT,OHDR - outputs
;
; OPTIONAL KEYWORDS:
;	soho - Boolean indicating SoHO (rather than Earth) view
;	date - A date specifier to override the header
;	t0 - If specified, causes the image(s) to be differentially derotated 
;		to match the specified time. 
;	lon - a range of longitudes to use
;	lat - a range of latitudes to use
;	lstep - The pixel size, in degrees of lon-/lat- itude, for the output
;		image at disk center.  If not specified, pixels at disk 
;		center are kept constant size in the transformation.  
;		If this is a 2-array, the (lon,lat) size is specified. 
;	quick - Causes all the images to be treated as happening at the same 
;		time, for P0, B0, and R0 calculation purposes.  Useful for 
;		faster processing of batches of images.
;	howard - Use the Howard et al rotation rate for small magnetic 
;		features (option to diff_rot) instead of the default Allen 
;		numbers.
;	rigid - Use the differential rotation speedf or the center of the 
;		rotated map, but treat the sun as a rigid body.  Boolean.
;	rlat - Specify the latitude whose rotation rate is to be used
;		for rotation.  Implies "rigid".
;	carrington - Boolean.  Use the Carrington rotation speed.  Implies
;		"rigid".
;
; HISTORY - Written 30-Apr-97, Craig DeForest
;	    Added `howard' keyword for pass-thru to zhelio2xy, 13-may-97 CED
;	    Added "Derotated-to:" to comment of output header, 22-may-97 CED
;	    Fixed default value of sine keyword, 3-Aug-97 CED
;	    Modified for proper CRPIXn usage, 13-Mar-98 CED
;	    Several minor replairs in vectorization etc., 15-Jun-98
;		(CED, Clare Parnell, and Daniel Brown)
;	    Added /rigid, 19-July-1998 CED
;	    Modified lstep to allow vector lstep specification, 16-Feb-99 CED
;	    Added VISIBLE in call to ZHELIO2XY, 16-Feb-99 CED
;	    Fixed auto-ranging of lat and lon, 17-Feb-99 CED
;	    Added ZGUESS_LOCATION call, 17-Feb-99 CED
;	    Miscellaneous fixes, plus missing-value setting is now allowed.
;	    Backed out the IDL-5 square-bracket stuff (some international
;	    colleagues are still using IDL-4)., 10-May-99, CED
;	    Fixed some vectorization stuff for auto-ranging, 9-Jun-99, CED
;	    Fixed $#(*@!! integer arithmetic in header calculation that was
;		causing half-pixel offsets in some maps. 27-Jul-99, CED
;
;-
pro ZHELIOGRAPHINIZE,in=in,out=out,ihdr=ihdr,ohdr=ohdr,lon=lon,lat=lat,lstep=lstep1,verbose=v,quick=quick,soho=soho,location=location,t0=t0,sine=sine,howard=howard,rigid=rigid,rlat=rlat,carrington=carrington,lerp=lerp,missing=missing

if not(isvalid(v)) then v=0
if not isvalid(missing) then missing= 0 * in(0,0,0)

if not isvalid(soho) then begin
	if not isvalid(location) then location = zguess_location(ihdr(0))
	if(strupcase(location) eq 'SOHO') then soho=1 else soho=0
end

sine = keyword_set(sine)

if(v) then print,"zheliographinize..."
if(v) then help,in

if (isvalid(t0)) then begin
	it0 = t0
	if (t0 eq fix(t0)) and (0 ge t0) and (t0 lt n_elements(ihdr)) then begin
		t0 = ihdr(t0).DATE_OBS + ' ' + ihdr(t0).TIME_OBS
		print,"ZHELIOGRAPHINIZE: Converted t0 index ",fug(it0)," to time-string '",t0,"'"
	end
end

i2hdr = ihdr
nim = n_elements(i2hdr)

zhdrunit,i2hdr,'solar-radii'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Auto-ranging:  if lon and lat ranges aren't specified, we need to 
; pick a set.  We generate two bogus maps -- one with latitude and one
; with longitude -- and zdeheliographinize them to the shape of the current
; image, then pick off the max and min.  Kludgy but it works.
;
; We do this every time now, because (A) it's not all that expensive, and
; (B) it lets us skip images in large pipelines with no intersection.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

if isvalid(t0) then t00 = t0

; Bogus maps are one degree spacing...
lons = (indgen(360) - 180) # (fltarr(180)+1)
lats = (fltarr(360) + 1)   # (indgen(180) - 90)

; Generate bogus-map header
lhdr = ihdr(0)
lhdr.naxis1 = 360
lhdr.naxis2 = 180
lhdr.crpix1 = 180 + 1
lhdr.crpix2 = 90 + 1
lhdr.crval1 = 0
lhdr.crval2 = 0
lhdr.ctype1 = 'degrees-longitude'
lhdr.ctype2 = 'degrees-latitude'
lhdr.cdelt1 = 1.0
lhdr.cdelt2 = 1.0
lhdr.date_obs = '1-Jan-1999'

; Generate destination header for zdeheliographinize; maximum of 10,000 points.
ilhdr = ihdr(0)
scale = sqrt(i2hdr(0).naxis1 * i2hdr(0).naxis2) / 30 ; Aim for a 30x30 final image

if(scale gt 1.0) then begin
	ilhdr.naxis1 = ilhdr.naxis1 / scale
	ilhdr.naxis2 = ilhdr.naxis2 / scale
	ilhdr.crpix1 = ilhdr.crpix1 / scale
	ilhdr.crpix2 = ilhdr.crpix2 / scale
	ilhdr.cdelt1 = ilhdr.cdelt1 * scale
	ilhdr.cdelt2 = ilhdr.cdelt2 * scale
end

zdeheliographinize,in=lons,ihdr=lhdr,out=lon2,ohdr=lon2hdr,tohdr=ilhdr
zdeheliographinize,in=lats,ihdr=lhdr,out=lat2,ohdr=lat2hdr,tohdr=ilhdr

; Check for valid lat/lon specs and fill in the blanks if necessary
lon0 = (-90)>([min(lon2),max(lon2)]+[-10,10])<(90)
lat0 = (-90)>([min(lat2),max(lat2)]+[-10,10])<(90)

if(n_elements(lon) ne 2) then lon=lon0
if(n_elements(lat) ne 2) then lat=lat0

; Check whether there are *any* valid points to convert!
if  (max(lon0) lt min(lon)) or (min(lon0) gt max(lon)) or $
    (max(lat0) lt min(lat)) or (min(lat0) gt max(lat)) then dummyflag = 1

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; lstep calculation:  If lstep isn't specified, we calculate it so that
; all our steps are the same as the pixels at disk center in the image.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if isvalid(lstep1) then lstep = lstep1 else lstep = [i2hdr(0).CDELT1,i2hdr(0).CDELT2] * (180/!PI)
if nlm(lstep) eq 1 then lstep = [lstep,lstep]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; rigid-body rotation...
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if keyword_set(rigid) and not isvalid(rlat) then rlat = (lat(0)+lat(1))/2.0
if keyword_set(carrington) or isvalid(rlat) or keyword_set(rigid) then rigid1 = 1
if isvalid(rlat) and keyword_set(carrington) then message,"You can't simultaneously set 'rlat' and '/carrington'."

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Generate the new header
;
ohdr = i2hdr

ohdr.NAXIS1 = (lon(1)-lon(0))/lstep(0)
ohdr.CROTA   = 0
ohdr.CRVAL1 = (lon(0)+lon(1))/2.0

if(sine) then begin
	ohdr.CTYPE1 = 'degrees-longitude'
	ohdr.CDELT1 = lstep(0)
	ohdr.CRPIX1 = ohdr.NAXIS1/2.0 + 1

	ohdr.NAXIS2 = !radeg * (sin(lat(1)/!radeg)-sin(lat(0)/!radeg))/float(lstep(1))
	ohdr.CRVAL2 = (sin(lat(0)/!radeg)+sin(lat(1)/!radeg))/2.0

	ohdr.CTYPE2='sine-latitude'
	ohdr.CDELT2 = lstep(1) / !radeg ; radeg converts to radians latitude!
	ohdr.CRPIX2 = ohdr.NAXIS2/2.0 + 1
end else begin
	ohdr.CTYPE1 = 'degrees-longitude'
	ohdr.CDELT1 = lstep(0)
	ohdr.CRPIX1 = ohdr.NAXIS1/2.0 + 1

	ohdr.NAXIS2 = (lat(1)-lat(0))/lstep(1)
	ohdr.CRVAL2 = (lat(0)+lat(1))/2.0
	ohdr.CTYPE2 = 'degrees-latitude'
	ohdr.CDELT2 = lstep(1)
	ohdr.CRPIX2 = ohdr.NAXIS2/2.0 + 1
end

if(nim gt 1) then begin
	out = replicate(in(0),ohdr(0).NAXIS1, ohdr(0).NAXIS2, nim)
	buf = replicate(in(0),ohdr(0).NAXIS1,ohdr(0).NAXIS2)
	for i=0,nim-1 do begin
		if(v gt 0) then print,"*** ZHELIOGRPHINIZE:  Recursing for image #",fug(i)," of ",fug(nim),"..."
		zheliographinize,in=in(*,*,i),ihdr=i2hdr(i),out=buf,ohdr=hbuf,lon=lon,lat=lat,lstep=lstep,verbose=v,q=q,t0=t0,sine=sine,soho=soho
		out(*,*,i) = buf
		ohdr(i) = hbuf

		if(v gt 0) then print,'* (Returned from recursion for image ',fug(i),' of ',fug(nim),')'
	end

	return
end
dummy = zcheck_hdr(i2hdr,/fix)

if isvalid(t0) then begin
	w=where(ohdr.COMMENT eq '')
	if(w(0) gt -1) then ohdr.COMMENT(w) = "Derotated-to: '"+t0+"'"
	if (v) then print,"Derotated to "+t0
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Generate the new image
if( ((size(out))(1) ne ohdr(0).NAXIS1) or ((size(out))(2) ne ohdr(0).NAXIS2) ) then begin
	if(v) then print,"Allocating output image..."
	out = replicate(missing,ohdr(0).NAXIS1,ohdr(0).NAXIS2,n_elements(in(0,0,*)))
end else $
	out(*,*) = missing

if(v gt 0) then begin
  print,'Longitude ranging from ',fug(lon(0)),' to ',fug(lon(1))
  print,'Latitude  ranging from ',fug(lat(0)),' to ',fug(lat(1))
  print,'Lstep is ',lstep,'; final im size is ',fug(ohdr.NAXIS1),'x',fug(ohdr.NAXIS2)
  if isvalid(t0) then print,'Differential-rotating to time:',t0,' (header time is ',ohdr(0).DATE_OBS+' '+ohdr(0).TIME_OBS,')'
end

if keyword_set(dummyflag) then begin
	if(v) then print,'Aborting rapidly - no intersection of latitude and longitude...'
	return
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Now, fill in the final image, chunk by chunk...
chunksize = long(100)
i2hdr = ihdr

lix = lindgen(ohdr(0).NAXIS1 * chunksize) mod ohdr(0).NAXIS1  ; Some nice indexed *integer*
liy = lindgen(ohdr(0).NAXIS1 * chunksize) / ohdr(0).NAXIS1    ; arrays (to index pixels out)

if(sine) then 					    $
	lyt = liy*lstep(1) + !radeg*sin(lat(0)/!radeg) $
else  						    $
	lyt = liy*lstep(1) + lat(0)		

if(v) then print,"min(lyt)=",min(lyt),"; max(lyt)=",max(lyt)

lx = (lix)*lstep(0) + lon(0)
ly = fltarr(ohdr(0).NAXIS1 * chunksize)

for iy=0,(ohdr(0).NAXIS2)/chunksize do begin
	
	if(sine) then begin
		ly = !radeg * asin((-89.99>(iy*chunksize*lstep(1)+lyt)<89.99)/!radeg)
		if(v) then print,"		min(ly)=",min(ly),"; max(ly)=",max(ly)
	end else $
		ly(*) = -89.99>((iy*chunksize*lstep(1)) + lyt)<89.99

	if(v gt 0) then print,'chunk '+fug(iy)+' of '+fug(ohdr(0).NAXIS2/chunksize)+'; ly=',fug(ly(0)),'...',format='(A,A,A,A,A,A,A,$)'
	if(v gt 1) and isvalid(date) then print,'; date=',date,format='(A,A,$)'
	if(v gt 1) then print,'; hdr.date_obs=',i2hdr.date_obs,'; i2hdr.time_obs=',i2hdr.time_obs,format='(A,A,A,A,$)'

	xy = zhelio2xy(ly,lx,i2hdr(0),b=b,l=l,p=p,date=date,soho=soho,verbose=(v gt 2),t0=t0,howard=howard,rigid=rigid,rlat=rlat,carrington=carrington,visible=visible)
	w = where(xy(0,*) ge 0 and xy(0,*) lt ihdr(0).NAXIS1 and $
	 	xy(1,*) ge 0 and xy(1,*) lt ihdr(0).NAXIS2 and visible)
		if(n_elements(w) gt 1) then $
			if(keyword_set(lerp)) then $
				out(lix(w),liy(w)+chunksize*iy) $
				 = interpolate(in,xy(0,w),xy(1,w)) $
			else $
				out(lix(w),liy(w)+chunksize*iy) $
				 = in(xy(0,w)+0.5,xy(1,w)+0.5) 
	delvarx,xy
	if(v gt 0) then begin
		print,' ',fug(n_elements(w)),' points'
	end

end

end