pro ZHOVER,in=in,ihdr=ihdr,out=out,ohdr=ohdr,latlon=latlon,scale=scale,naxis=naxis,t0=t0,howard=howard,rigid=rigid,rlat=rlat,lerp=lerp,verbose=v,missing=missing

if not isvalid(missing) then missing = 0.0

if keyword_set(rigid) and not isvalid(rlat) then rlat = latlon(0)

thdr = ihdr(0)
thdr.naxis1 = naxis
thdr.naxis2 = naxis
thdr.crpix1 = naxis/2.0
thdr.crpix2 = naxis/2.0
thdr.crval1 = 0
thdr.crval2 = 0
thdr.ctype1 = 'solar-radii'
thdr.ctype2 = 'solar-radii'
thdr.cdelt1 = scale
thdr.cdelt2 = scale
thdr.HISTORY[0]="NOTE REPROJECTED B0,L!"
thdr = add_tag(temporary(thdr),latlon(0),'B0')
thdr.B0 = latlon(0)
thdr = add_tag(temporary(thdr),latlon(1),'L')
thdr.L = latlon(1)
thdr.time_obs = ''
thdr.date_obs = t0

out = fltarr(naxis,naxis,/nozero)
out(*) = missing

;;
;; Calculate lat/lon extend by heliographinizing the *output* image plane, to 
;; a coarse resolution  This lets zheliographinize determine the lat/lon extent
;; that are required...
;;

print,'Calculating lat/lon extent...'
foo = replicate(1b,thdr.naxis1,thdr.naxis2)
zheliographinize,in=foo,ihdr=thdr,ohdr=fhdr,t0=t0,howard=howard,rigid=rigid,rlat=rlat,lstep=1.0,missing=0

slop = 4

lonrange = ((fhdr.crval1 + ([0,fhdr.naxis1] - fhdr.crpix1) * fhdr.cdelt1 + [-1.0,1.0]*slop  + 180) + 7200) mod 360 - 180
latrange = fhdr.crval2 + ([0,fhdr.naxis2] - fhdr.crpix2) * fhdr.cdelt2 + [-1,1] * slop  < 90 > (-90)

;lonrange=latlon(1)+[-20,20]
;latrange=latlon(0)+[-20,20]
print,'lonrange=',lonrange, 'latrange=',latrange


for i=0,nlm(ihdr)-1 do begin
    print,"Handling layer #",i
    zheliographinize,in=in(*,*,i),ihdr=ihdr(i),out=im0,ohdr=hdr0,t0=t0,howard=howard,rigid=rigid,rlat=rlat,lon=lonrange,lat=latrange,lerp=lerp,v=v,missing=missing
print,"    zdehelio..."
    zdeheliographinize,in=im0,ihdr=hdr0,b=latlon(0),l1=latlon(1),tohdr=thdr,out=im1,ohdr=hdr1,lerp=lerp,v=v,missing=missing
    if(i eq 0) then begin
        out =  im1
        ohdr = hdr1
    end else begin
        out = [[[out]],[[im1]]]
        ohdr = [ohdr,hdr1]
    end
end
end