;+
; PURPOSE:
;  cartesian to log pol transform
;      We take the log of the elongation to be able to display all the
;      instruments on the same image.
;
; INPUTS:
;  h : array of pointers to wcs compliant headers. Note that the images will 
;      be layered in that order
;  im : array of pointers to images, in the same order than for h
;  imsizeout : [x,y] size of the output image, x axis is log elong, y axis is pa
;  elongminmax : [min,max] elongation in degrees
;  paminmax : [min,max] heliographic polar angle in degrees (0: north pole, 90: east limb, 180: south pole, 270: west limb)
;  imelonminmax : array of pointers to 2 elements vectors. The pointer order should be the same than for h. 
;                 Each 2 elements vector contain the min and max elongation of each image to be transformed. (deg)
;  linear : if set use a linear scale for elong. 
;
; OUTPUTS:
;  The function returns an image.
;  wcsout : wcs structure describing the image. Note that the SLG projection is not a wcs standard. We should propose that it becomes part of the standard someday.
;
; EXAMPLE:
;  impol=scccart2logpol(im,h,[1024,600],[1.5,88],[0.,120])
;
; See also W.T.Thompson, A&A Dec 15, 2005, section 4 for more details about the helioprojective coordinates.
;-

function scccart2logpol,im,h,imsizeout,elongminmax,paminmax,wcsout=wcsout,imelongminmax=imelongminmax,linear=linear

; ---- compute coordinates for each of the transformed image
if keyword_set(linear) then elong=replicate(1,imsizeout[1])##lgen(imsizeout[0],elongminmax) else elong=replicate(1,imsizeout[1])##(10^lgen(imsizeout[0],alog10(elongminmax)))
pa=lgen(imsizeout[1],paminmax)##replicate(1,imsizeout[0])

; ---- init output image
imo=make_array(dim=imsizeout,/float)

; ---- convert in cartesian heliocentric coords
polcoord=fltarr(2,imsizeout[0],imsizeout[1])
polcoord[0,*,*]=elong
polcoord[1,*,*]=pa
cartcoord=heliorad2heliocart(reverse(polcoord,1))

; ---- test if input is a single image or a list of images
; if size(h,/type) ne 10 then begin
; 	; -- seems like it's a single image	
; endif	

; ---- init default elongation range 
elongrange=elongminmax


nbim=n_elements(h)
; ---- compute transform for each image
for i=0,nbim-1 do begin
	; ---- get image heliocentric cart pixel position
	if tag_exist(*h[i],'coord_type') then wcs=*(h[i]) else wcs=fitshead2wcs(*(h[i]),system='B')
	; -- use proper units
	;fw2deg=(wcschangeunits(wcs.cunit[0],'deg',1.))[0]
	fdeg2w=(wcschangeunits('deg',wcs.cunit[0],1.))[0]
	
	; ---- get the pixel position 
	pixpos=wcs_get_pixel(wcs,cartcoord*fdeg2w,/force_proj)
	
	; ---- interpolate
	imotmp=interpolate(*(im[i]), reform(pixpos[0,*,*]), reform(pixpos[1,*,*]),cubic=-0.5,missing=-32000.)

	if n_elements(imelongminmax) gt 0 then elongrange=*(imelongminmax[i])

	mok=where((imotmp ne -32000.) and (elong ge elongrange[0]) and (elong le elongrange[1]),nbok)
	
	if nbok gt 0 then imo[mok]=imotmp[mok]
	
endfor

; ---- create a header describing the ouput image
if keyword_set(linear) then cdelt1=(elongminmax[1]-elongminmax[0])/(imsizeout[0]-1) else cdelt1=(alog10(elongminmax[1])-alog10(elongminmax[0]))/(imsizeout[0]-1)
cdelt2=(paminmax[1]-paminmax[0])/(imsizeout[1]-1)
if keyword_set(linear) then crpix=[1.,1.] else crpix=[-alog10(elongminmax[0])/cdelt1,1.] ; !!! fortran convention !!!

wcsout={coord_type:'Helioprojective-Radial',$
	wcsname:'Helioprojective-Radial',$
	naxis:imsizeout,$
	variation:'PC',$
	compliant:1,$
	projection:'SLG',$
	ix:0L,$
	iy:1L,$
	crpix:crpix,$  ; !!! fortran convention !!!
	crval:[polcoord[1,0,0],polcoord[0,0,0]-90.],$
	ctype:['SOLI-SLG','HCPA-SLG'],$
	cname:['',''],$
	cunit:['deg','deg'],$
	cdelt:[cdelt1,cdelt2],$
	pc:double(identity(2)),$
	proj_names:['LONPOLE','PV2_1','PV1_1','P1_2'],$
	proj_values:[180d,0,polcoord[1,0,0],polcoord[0,0,0]-90.],$
	roll_angle:0d,$
	simple:1B,$
	time:0,$
	position:0}

return,imo
end