;+ ; NAME: ; render ; PURPOSE: ; Render a rectangular projection map to a sphere. ; DESCRIPTION: ; The most confusing aspect of this program is by far the image orientation. ; I can't make the conventions natural for everyone (not even myself), but ; here's what this program does: ; ; In the output image there is a natural coordinate system, (x,y). Each ; row, ie., image[*,i] corresponds to a fixed value of y. image[*,0] ; is the row with the LOWEST value of y. This may seem strange, but it ; permits a "normal" orientation of the image if displayed with !order=0 ; (the usual IDL default). The x direction (image[i,*]) works in the ; obvious way, image[0,*] is the column with the LOWEST value of x. So, ; viewed with !ORDER=0, this is an image that has up at the top. The ; top is also considered NORTH when dealing with images that one sees in ; the sky (or relative to your point of view). ; ; The position angle of the pole sets the rotation of the sphere in the ; plane of the sky (image). North is to the top (+y) and East is to the ; left (-x). The position angle (POLE) is the angle of the apparent ; North pole of the sphere relative to up (North) in the image, measured ; From up (North) counter-clockwise (also known as eastward from North). ; ; Visualizing the indexing of the input map is equally confusing. In ; the coordinate system containing the map, North is +y, East is -x, as ; before. These coordinates are plane-of-sky on the object. The map ; has its own coordinates, latitude and longitude. The first index into ; the map is treated like x but is actually East Longitude. The left edge ; of image[0,*] is precisely at 0 degrees longitude. As the first index ; increases (increasing x?), the longitude increases. The right edge of ; the last column in the map is 360 degrees. ; ; Ok, here is the weird part. The first row in the map is the NORTHERN-most ; set of pixels. The last row is the SOUTHERN-most set. So, map[0,0] ; touches the prime meridian and the north pole. ; Thus map[i,j]: ; longitude = (i+0.5)/nl * 360.0 ; latitude = 90.0 - (j+0.5)/nt * 180.0 ; where nl is the width of map and nt is the height of the map. ; These formulas give the lat,lon of the CENTER of the pixel in degrees. ; ; Now, understanding the sub-"earth" and sub-solar longitudes becomes ; easy. The sub-"earth" point is the lat,lon of the map that is in the ; center of the projected disk. The sub-solar point is the lat,lon of ; the map that is nearest to the sun (normal solar illumination). Note, ; that as the object rotates, the sub-earth longitude will decrease, not ; increase. ; ; CATEGORY: ; Image display ; CALLING SEQUENCE: ; render,map,radius,scale,pole,selat,selon,sslat,sslon,nx,ny,image ; INPUTS: ; map - Array containing a full map of surface ; radius - Radius of object (eg., in km) ; scale - Scale of image (eg., km/pixel) ; pole - Position angle of pole, east from north (degrees) ; selat - Sub-earth latitude (degrees) ; selon - Sub-earth longitude (degrees) ; sslat - Sub-solar latitude (degrees). Only used for Hapke functions. ; sslon - Sub-solar longitude (degrees). Only used for Hapke functions. ; nx - X size of output image in pixels. ; ny - Y xize of output image in pixels. ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; SILENT - Flag, if set suppresses all printed output ; NODISPLAY - Flag, if set suppresses graphical output ; GEOM - Undefined (default), no special action. ; If defined but not a structure, then upon return geom will contain ; all the geometric information needed to render the image with a ; new map. ; If defined and is a structure, then the contents are taken to be ; that needed to do the final image calculation without redoing ; all the geometric computation. ; ; Limb darkening model to use (set one of these, at most): ; ; LIN - Linear limb-darkening coeff, map is normal albedo ; MIN - Minnaert limb-darkening coeff, map is normal albedo ; HAPKE - Hapke scattering model, map is single scattern albedo ; HAPKE[0] = h (old style, circa 1981) ; HAPKE[1] = P(0) ; HAP2 - 2nd generation Hapke scattering model, 1986 and the book, ; "Theory of Reflectance and Emittance Spectroscopy" ; HAP2[0] = h (new style, circa 1986, p. 226 in book) ; HAP2[1] = P(g) ; HAP2[2] = B0, Emperical backscatter factor (p. 228) ; HAP2[3] = Theta(bar), surface roughness parameter. ; H93 - When used with HAP2, will select which approximation to the ; H-function is used. Default is the simplest (and fastest) ; formula from the 1981 Hapke paper. H93, when set, will ; force the use of eqn. 8.57 from the Hapke book (from 1993). ; ; OUTPUTS: ; image - Rendered image of object. ; KEYWORD OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; 95/07/24, Written by Marc W. Buie, Lowell Observatory ; 96/04/11, MWB, fixed numerous bugs, now works for sun and viewpoint ; being different ; 97/08/21, MWB, added HAP2 keyword ; 97/08/27, MWB, added GEOM keyword ; 97/09/18, MWB, repaired incorrect mapping of map to surface. ;- pro render,map,radius,scale,in_pole,selat,selon,sslat,sslon,nx,ny,image, $ LIN=xlimb,MIN=k,HAPKE=hapke,HAP2=hapke2,NODISPLAY=nodisplay, $ SILENT=silent,GEOM=geom,H93=h93,XARR=xarr,YARR=yarr ; Save size of input map sz=size(map) nlon = sz[1] nlat = sz[2] ; Display input map if not keyword_set(nodisplay) then begin setwin,1,xsize=nlon,ysize=nlat tvscl,map,/order endif ; Set up output image array image=fltarr(nx,ny) ; Check on the GEOM keyword and find out what to do. szg=size(geom) typecode = szg[n_elements(szg)-2] fullcompute = n_elements(geom) le 1 IF fullcompute THEN BEGIN ; Set up image that contains the y coordinate of each pixel (relative to cntr) y = ((findgen(ny) - float(ny-1)/2.0) * scale)/radius onex = replicate(1.0,nx) yarr = onex#y ; Set up image that contains the x coordinate of each pixel (relative to cntr) x = ((findgen(nx) - float(nx-1)/2.0) * scale)/radius oney = replicate(1.0,ny) xarr = x#oney ; Image of pixel distance from center. rsq = xarr^2 + yarr^2 ; Convert input angles to radians pole = in_pole / !radeg ; position angle of pole theta = selat / !radeg ; sub-earth latitude phi = selon / !radeg ; sub-earth longitude stheta = sslat / !radeg ; sub-solar latitude sphi = sslon / !radeg ; sub-solar longitude ; Unit vector that points toward observer xv = cos(theta)*cos(phi) yv = cos(theta)*sin(phi) zv = sin(theta) ; Unit vector that points toward sun xs = cos(stheta)*cos(sphi) ys = cos(stheta)*sin(sphi) zs = sin(stheta) ; Index into IMAGE that fall on disk of object. zon = where(rsq lt 1.0) ; 3-D coordinates of image plane ON the object. (this is also the set ; of unit vectors for each of the points on the surface of the object) x = xarr[zon] y = yarr[zon] z = sqrt( 1.0 - rsq[zon] ) ; Rotate the surface vectors to the viewpoint rotpoint,x,y,z,'z',-pole,xp,yp,zp rotpoint,xp,yp,zp,'x',-theta,xpp,ypp,zpp rotpoint,xpp,ypp,zpp,'y',phi,xppp,yppp,zppp ; Compute dot products with viewpoint and illumination directions mu = zppp*xv + xppp*yv + yppp*zv smu = zppp*xs + xppp*ys + yppp*zs g = acos(xv*xs + yv*ys + zv*zs) ; Compute map coordinates for each surface pixel colat = acos(yppp) lon = atan(xppp,zppp) ; Convert negative longitudes to positive. zl = where(lon lt 0,count) if count ne 0 then lon[zl]=lon[zl]+!pi*2.0 ; Convert from lat,lon to indicies into the map, these arrays now point ; to the visbile pixels, but they are not necessarily illuminated. xidx = long(lon/!pi/2.0*nlon) yidx = long(colat/!pi*nlat) ; colat can be exactly !pi, must map that back down one pixel zc = where(yidx eq nlat, count) if count ne 0 then yidx[zc]=nlat-1 ; Find only those pixels that are illuminated ill = where(smu gt 0.,ill_count) ; compute image pixels for relevant reflectance "law" if ill_count gt 0 then begin ; Linear limb darkening case (Marcialis) if keyword_set(xlimb) then begin if not keyword_set(silent) then print,'Linear limb darkening' image[zon[ill]] = map[xidx[ill],yidx[ill]]*(1.0-xlimb+xlimb*mu[ill]) ; Minnaert limb darkening endif else if keyword_set(k) then begin if not keyword_set(silent) then print,'Minnaert limb darkening' image[zon[ill]] = map[xidx[ill],yidx[ill]]*mu[ill]^(2*k-1.0) ; Hapke scattering function. endif else if keyword_set(hapke) then begin if not keyword_set(silent) then print,'Hapke theory used.' image[zon[ill]] = bidr(map[xidx[ill],yidx[ill]],mu[ill],smu[ill],g[ill], $ hapke[0],hapke[1])*!pi ; Hapke scattering function, second generation. endif else if keyword_set(hapke2) then begin if not keyword_set(silent) then print,'Hapke theory used, new formalism.' image[zon[ill]] = bidr2(map[xidx[ill],yidx[ill]],mu[ill],smu[ill],g[ill], $ hapke2[0],hapke2[1],hapke2[2],hapke2[3],h93=h93)*!pi ; No limb darkening at all, just wrap the map as is. endif else begin if not keyword_set(silent) then print,'No limb darkening used.' image[zon[ill]] = map[xidx[ill],yidx[ill]] endelse IF typecode ne 0 THEN BEGIN geom = [xidx[ill], $ yidx[ill], $ mu[ill], $ smu[ill], $ g[ill], $ zon[ill] ] ENDIF endif ENDIF ELSE BEGIN ; GEOM has geometry illc = n_elements(geom)/6 xidx = 0 yidx = 1 mu = 2 smu = 3 g = 4 zon = 5 ; compute image pixels for relevant reflectance "law" ; Linear limb darkening case (Marcialis) if keyword_set(xlimb) then begin if not keyword_set(silent) then print,'Linear limb darkening' image[geom[zon*illc:(zon+1)*illc-1]] = $ map[geom[xidx*illc:(xidx+1)*illc-1], $ geom[yidx*illc:(yidx+1)*illc-1]] * $ (1.0-xlimb+xlimb*geom[mu*illc:(mu+1)*illc-1]) ; Minnaert limb darkening endif else if keyword_set(k) then begin if not keyword_set(silent) then print,'Minnaert limb darkening' image[geom[zon*illc:(zon+1)*illc-1]] = $ map[geom[xidx*illc:(xidx+1)*illc-1], $ geom[yidx*illc:(yidx+1)*illc-1]] * $ geom[mu*illc:(mu+1)*illc-1]^(2*k-1.0) ; Hapke scattering function. endif else if keyword_set(hapke) then begin if not keyword_set(silent) then print,'Hapke theory used.' image[geom[zon*illc:(zon+1)*illc-1]] = $ bidr(map[geom[xidx*illc:(xidx+1)*illc-1], $ geom[yidx*illc:(yidx+1)*illc-1]], $ geom[mu*illc:(mu+1)*illc-1], $ geom[smu*illc:(smu+1)*illc-1], $ geom[g*illc:(g+1)*illc-1], $ hapke[0],hapke[1])*!pi ; Hapke scattering function, second generation. endif else if keyword_set(hapke2) then begin if not keyword_set(silent) then print,'Hapke theory used, new formalism.' image[geom[zon*illc:(zon+1)*illc-1]] = $ bidr2(map[geom[xidx*illc:(xidx+1)*illc-1], $ geom[yidx*illc:(yidx+1)*illc-1]], $ geom[mu*illc:(mu+1)*illc-1], $ geom[smu*illc:(smu+1)*illc-1], $ geom[g*illc:(g+1)*illc-1], $ hapke2[0],hapke2[1],hapke2[2],hapke2[3],h93=h93)*!pi ; No limb darkening at all, just wrap the map as is. endif else begin if not keyword_set(silent) then print,'No limb darkening used.' image[geom[zon*illc:(zon+1)*illc-1]] = $ map[geom[xidx*illc:(xidx+1)*illc-1], $ geom[yidx*illc:(yidx+1)*illc-1]] endelse ENDELSE ; Display final image. if not keyword_set(nodisplay) then begin setwin,0,xsize=nx,ysize=ny tvscl,image endif end