;+ ; NAME: ; PRO VISU_CUBE ; ; PURPOSE: ; Create a cube of images for a series of epochs in the fits format ; ; CALLING SEQUENCE: ; VISU_CUBE, [workDirectory=workDirectory] [, targetName=targetName] [, img_reso=img_reso] [, viewfield=viewfield] ; [, nbd=nbd] [, epoch=epoch] [, observer=observer] [, diameter=diameter] [, magnification=magnification] ; [, outfile=outfile] [, timescale=timescale] [, typeCoord=typeCoord] ; ; INPUT: ; ; ; ; KEYWORD PARAMETERS: ; WORKDIRECTORY - string specifying the work directory where files will be stored (default : '') ; TARGETNAME - string name of the object to be visualized (default : lutetia) ; IMG_RESO - dimensions of the images to be plotted [int X_Dim, int Y_Dim] (default : [2048, 2048]) ; VIEWFIELD - dimensions of the field of view in degrees. If not specified, the field of view is adapted to the dimensions of the targeted object ; NBD - number of images to be plotted (default : 1) ; EPOCH - epoch of the first view ; OBSERVER - string name of the space probe (default : rosetta) ; DIAMETER - diameter of the imaging instrument for psf computation in meters (default : 2.0) ; MAGNIFICATION - magnification of the imaging instrument for psf computation (default : 1.0) ; OUTFILE - string name of the file to be plotted (should end by '.fits') ; TIMESCALE - default : 'UTC' ; TYPECOORD - default : 'Apparent of date' ; ; OUTPUT: ; fits data cube ; ; CREATION ; 2012-mar-08: J-B.Valet (IMCCE, Observatoire de Paris) ; ; ; HISTORY ; ;- pro visu_cube, $ workDirectory=workDirectory, $ targetName=targetName, $ targetNumber=targetNumber, $ img_reso=img_reso, $ ;1x2 array - dimensions of image (nb of pixels along the x and y axis) viewfield=viewfield, $ ;field of view of the instrument (DEG) nbd=nbd, $ ;nb of dates to be computed step=step, $ epoch=epoch, $ ;starting epoch (julian day) observer=observer, $ ;string outfile=outfile, $ ;string mapFile=mapFile, $ ;string pointing an image file psf=psf, $ ;string pointing a fits file airy=airy, $ ;boolean seeing=seeing, $ ;float aperture=aperture, $ ;float wavelength=wavelength, $ ;float shademod=shademod, $ ;string ('polyshade' or 'lsl' for a more complex computation of light diffusion) opt=opt, $ RV=RV, $ alt=alt, $ stars=stars, $ fconfig=fconfig COMMON share, objVertex, objPolygons, barr, nbFacet if not keyword_set(workDirectory) then workDirectory='' ;data about the target if not keyword_set(targetName) then targetName='lutetia' if not keyword_set(targetNumber) then targetNumber='21' ;definition of the observer and used instrument if not keyword_set(observer) then observer='rosetta' if not keyword_set(img_reso) then img_reso=[512, 512] if not keyword_set(viewfield) then viewfield=[2.35, 2.35] ;default queried epochs if not keyword_set(nbd) then nbd=150 if not keyword_set(epoch) then epoch=2455388.1D if not keyword_set(step) then step=0.0007D if not keyword_set(outfile) then outfile='lut.fits' ; if not keyword_set(mapFile) then mapFile='mar1kuu2.tif' ;psf used ; if not keyword_set(seeing) then seeing = 50. ; if not keyword_set(psf) then psf = 'psf.Ks.01037-reduc.fits' if not keyword_set(psf) then airy=1 if not keyword_set(aperture) then aperture = 0.01D if not keyword_set(wavelength) then wavelength = 1000.E-8 ;light diffusion model if not keyword_set(shademod) then shademod = 'polyshade' ; if not keyword_set(shademod) then shademod = 'lsl' if not keyword_set(stars) then stars = 1 ;visualizations if not keyword_set(fconfig) then begin if not keyword_set(RV) then RV = 0 if not keyword_set(opt) then opt = 1 if not keyword_set(alt) then alt = 0 endif vizierRadius = sqrt(TOTAL(viewfield*viewfield)) / 2 nbFrames = 0 if alt eq 1 then nbFrames++ if RV eq 1 then nbFrames++ if opt eq 1 then nbFrames++ nbFrames = nbd * nbFrames ; ; for i=0, n_elements(psf_head) - 1 do begin ; ; line = psf_head[i] ; ; if (strpos(line, 'CD1_1') ne -1) then begin ; ; line = strmid(line, 10) ; ; line = strtrim(strmid(line, 0, strpos(line, '/')), 2) ; ; PSF_CD11 = float(line) ; ; endif else if (strpos(line, 'CD1_2') ne -1) then begin ; ; line = strmid(line, 10) ; ; line = strtrim(strmid(line, 0, strpos(line, '/')), 2) ; ; PSF_CD12 = float(line) ; ; endif else if (strpos(line, 'CD2_1') ne -1) then begin ; ; line = strmid(line, 10) ; ; line = strtrim(strmid(line, 0, strpos(line, '/')), 2) ; ; PSF_CD21 = float(line) ; ; endif else if (strpos(line, 'CD2_2') ne -1) then begin ; ; line = strmid(line, 10) ; ; line = strtrim(strmid(line, 0, strpos(line, '/')), 2) ; ; PSF_CD22 = float(line) ; ; endif ; ; endfor ; ; PSF_CD = [[PSF_CD11, PSF_CD12], [PSF_CD21, PSF_CD22]] IMG_CD = [[viewfield[0] / img_reso[0], 0], [0, viewfield[1] / img_reso[1]]] ; if ((ABS(PSF_CD[0, 0] - IMG_CD[0, 0]) ne 0.) or (ABS(PSF_CD[1, 1] - IMG_CD[1, 1]) ne 0.)) then psf = psf_resample(psf, PSF_CD, IMG_CD) tmat = [IMG_CD[0, 0], IMG_CD[1, 1]] if keyword_set(psf) then begin print, 'psf user' psf = readfits(psf, psf_head) psf = psf / total(psf) endif else if keyword_set(airy) then begin print, 'psf airy' psf = psf_airy(tmat, aperture=aperture, wavelength=wavelength) psf = psf / total(psf) endif else if keyword_set(seeing) then begin print, 'psf seeing' psf = psf_seeing(tmat, seeing) psf = psf / total(psf) endif degToMas=3600000 ; frames = bytarr(img_reso(0), img_reso(1), nbFrames) ;images' cube (X_resolution x Y_resolution x Number_of_images) frames = fltarr(img_reso[0], img_reso[1], nbFrames) verfile='/astrodata/SHAPE_MODELS/VER/' + targetName + '.ver' ;----------------------------- Reading ver files and applying maps ------------------------------; ; for i=0, n_elements(verfile) do begin read_ver_shape_model, verFile, $ objVertex, objPolygons, $ nbVertex=nbVertex, $ nbFacet=nbFacet ; endfor albedoMap = fltarr(nbFacet) + 1. ;- FacetINFO( 0,*) : Number of the Facet First Vertex ;- FacetINFO( 1,*) : Number of the Facet Second Vertex ;- FacetINFO( 2,*) : Number of the Facet Third Vertex ;- FacetINFO( 3,*) : Facet Visible (1) or not (0) ;- FacetINFO( 4,*) : Reflecion angle ;- FacetINFO( 5,*) : Incidence angle ;- FacetINFO( 6,*) : Radial velocity toward the Earth ;- FacetINFO( 7,*) : Facet superficie ;- FacetINFO( 8,*) : Radial velocity toward the Sun ;- FacetINFO( 9,*) : Diffusion Law Value: Lambert ;- FacetINFO(10,*) : Diffusion Law Value: Lommel-Seelinger ;- FacetINFO(11,*) : Diffusion Law Value: Minnaert ;- FacetINFO(12,*) : Total intensity of the Facet ;- FacetINFO(13,*) : Albedo information (if any) altitudeMap = fltarr(nbFacet) + 1. if keyword_set(mapFile) then begin imageMap = READ_TIFF(mapFile) s = size(imageMap) if (s(1) eq 3) then imageMap=reform(imageMap[0, *, *]) s = size(imageMap) print, 'map facet' map = map_facets(objVertex, objPolygons, SPACING=(s(1) / 360), /EQUATORIAL) for k=0, nbFacet-1 do begin ;print, k area = imageMap(where(map eq k)) if (n_elements(area) ne 0) then begin m = total(area) / n_elements(area) albedoMap[k] = m endif endfor albedoMap = albedoMap * 255. / max(albedoMap) endif if keyword_set(alt) then begin altitudeMap = visu_altitude() endif ;------------------------------------- http request ephemph -------------------------------------; ;wget = 'wget -q -O ' + workfile + 'phEphemeris.vot' ; ;url = 'http://vo.imcce.fr/webservices/miriade/ephemph_query.php?' ; ;args = ['-type=aster', $ ; ;'-name=lutetia', $ ; '-name=' + name, $ ; '-ep=' + TRIM(STRING(epoch)), $ ; ;'-nbd=500', $ ; '-nbd=' + TRIM(STRING(nbd)), $ ; '-step=40.0s', $ ; '-tscale=UTC', $ ; '-observer=' + observer, $ ; ;'-teph=4', $ ; ;'-tcoor=1', $ ; ;'-rplane=1', $ ; '-mime=votable', $ ; '-view=none', $ ; '-from=MyProject'] ; ;opt = strtrim(args(0), 2) ; ;size_args = size(args) ; ;for i=1, size_args(1)-1 do begin ; opt = strtrim(opt, 2) + '&' + strtrim(args(i), 2) ;endfor ; ;cmd = wget + ' "' + strtrim(url, 2) + strtrim(opt, 2) + '"' ; ; ;print, 'envoi requete http' ; ;SPAWN, cmd, result ; ;--------------------------------- recuperation of ephemph result ---------------------------------; ; ;print, 'parsing du fichier d ephemerides physiques' ; ;parser = OBJ_NEW('phEphemerisParser') ;votFile = FILEPATH('phEphemeris.vot', $ ; SUBDIRECTORY = ['..', '..', '..', '..', '..', 'obs', 'jvalet']) ;parser->ParseFile, votFile ;sepLong = parser->GetSEPLong() ;sepLat = parser->GetSEPLat() ;sspLong = parser->GetSSPLong() ;sspLat = parser->GetSSPLat() ;npole = parser->GetNp() ;distG = parser->GetDg() ;distH = parser->GetDh() ;OBJ_DESTROY, parser ;------------------------------------- local request ephemph --------------------------------------; ;-- call to miriade ephemph cmd = 'ephemph aster -sonde Spice:rosetta -d ' + strtrim(string(nbd), 2) + ' -s 1 -print -p 0.0007 -j ' + strtrim(string(epoch), 2) + ' -nom ' + string(targetNumber) + ' -f ephem_ros_lut.vot' ; ; SPAWN, cmd, result, error ;-- open the ephemeris file from miriade ephemph openr, unit, 'Eproc/ephem/ephem_ros_lut.vot', /GET_LUN result='' i = 0 crval1=fltarr(nbd) crval2=fltarr(nbd) list = [] ;-- read the ephemph response in the file while ~ EOF(unit) do begin readf, unit, result if(strpos(result, '#') eq -1) and (strtrim(result) ne '') and (i lt nbd) then begin print, i ;-- Sub-solar and sub-Earth longitudes and latitudes sepLong = float(strmid(result, 22, 7)) ;SEP long sepLat = float(strmid(result, 29, 7)) ;SEP lat sspLong = float(strmid(result, 36, 7)) ;SSP long sspLat = float(strmid(result, 43, 7)) ;SSP lat ;-- Pole angle npole = float(strmid(result, 50, 7)) ;Np ;-- Distance of the target from the observer (distG) and from the Sun (distH) distG = float(strmid(result, 89, 13)) ;Dg distH = float(strmid(result, 102, 13)) ;Dh ;-- Right Ascension and declination of target as viewed from the observer RA = strtrim(strmid(result, 131, 15), 2) ;RA DE = strtrim(strmid(result, 146, 13), 2) ;DE ;-- convert the coordinates in degrees and store them in an array to be written in the fits header (CRVAL1, CRVAL2) res = RA crval1(i) = 0. p = strpos(res, ' ') crval1(i) = crval1(i) + strmid(res, 0, p) * 15.0 res = strtrim(strmid(res, p), 1) p = strpos(res, ' ') crval1(i) = crval1(i) + strmid(res, 0, p) * 0.25 crval1(i) = crval1(i) + strtrim(strmid(res, p), 1) / 240.0 res = DE crval2(i) = 0. if ((strpos(res, '+')) ne -1) then begin p = strpos(res, ' ') crval2(i) = crval2(i) + strmid(res, 0, p) * 1.0 res = strtrim(strmid(res, p), 1) p = strpos(res, ' ') crval2(i) = crval2(i) + strmid(res, 0, p) / 60.0 crval2(i) = crval2(i) + strtrim(strmid(res, p), 1) / 3600.0 endif else if ((strpos(res, '-')) ne -1) then begin crval2(i) = crval2(i) - strmid(res, 0, strpos(res, ' ')) * 1.0 res = strtrim(strmid(res, strpos(res, ' ')), 1) crval2(i) = crval2(i) - strmid(res, 0, strpos(res, ' ')) / 60.0 crval2(i) = crval2(i) - strtrim(strmid(res, strpos(res, ' ')), 1) / 3600.0 endif epoch = epoch + step * i if stars eq 1 then begin degCoord_cur = sexa2deg(RA, DE) ;-- If list don't exist --> initializing if n_elements(list) eq 0 then begin list = vizier(RA, DE, 1.4*vizierRadius) ;list of sources from the ppmx catalog (vizier request) degCoord_ori = degCoord_cur endif else begin d = sqrt((degCoord_ori[0] - degCoord_cur[0])^2 + (degCoord_ori[1] - degCoord_cur[1])^2) ;-- if the previous list isn't sufficient for the new field of view --> new vizier call if (d gt (0.2 * vizierRadius)) then begin list = vizier(RA, DE, 1.4 * vizierRadius) degCoord_ori = degCoord_cur endif endelse endif s = n_elements(list) backgrd = fltarr(img_reso[0], img_reso[1]) counter = 0 if opt eq 1 then begin visu, lsep=sepLong, $ bsep=sepLat, $ lsun=sspLong, $ bsun=sspLat, $ NP=npole, $ Dg=distG, $ Dh=distH, $ typecoord='Apparent of date', $ timescale='UTC', $ resolution=img_reso, $ verfile=verfile, $ aperture=2.0, $ viewField=viewfield, $ omega0=2., $ epoch=epoch, $ objVertex=objVertex, $ objPolygons=objPolygons, $ albedoMap=albedoMap, $ shademod=shademod, $ ; onlyVR=1, $ image=image if (stars eq 1) && (s gt 1) && (counter eq 0) then begin backgrd = draw_stars(RA, DE, viewfield, image, LIST=list) endif image = image + backgrd if keyword_set(psf) then image = psf_convolve(image, psf=psf) ; image = bytscl(image, MAX=1.0, MIN=0.0) frames[*, *, i] = image ;image = draw_label(image, img_reso, sepLong, sepLat, sspLong, sspLat, npole, targetName, epoch=epoch+0.0007*i, timescale) counter++ endif if RV eq 1 then begin print, 'visu_rv' visu_RV, lsep=sepLong, $ bsep=sepLat, $ lsun=sspLong, $ bsun=sspLat, $ NP=npole, $ Dg=distG, $ Dh=distH, $ typecoord='Apparent of date', $ timescale='UTC', $ resolution=img_reso, $ verfile=verfile, $ aperture=2.0, $ ;viewField=viewfield, $ omega0=2., $ epoch=epoch, $ objVertex=objVertex, $ objPolygons=objPolygons, $ albedoMap=albedoMap, $ shademod=shade;image = draw_label(image, img_reso, sepLong, sepLat, sspLong, sspLat, npole, targetName, epoch=epoch+0.0007*i, timescale)mod, $ ; onlyVR=1, $ image=image if (stars eq 1) && (s gt 1) && (counter eq 0) then begin backgrd = draw_stars(RA, DE, viewfield, image, LIST=list) endif image = image + backgrd if keyword_set(psf) then image = psf_convolve(image, psf=psf) image = bytscl(image, MAX=1.0, MIN=0.0) frames[*, *, i + counter] = image ;image = draw_label(image, img_reso, sepLong, sepLat, sspLong, sspLat, npole, targetName, epoch=epoch+0.0007*i, timescale) counter++ endif if alt eq 1 then begin print, 'visu_altitude' visu_altitude, lsep=sepLong, $ bsep=sepLat, $ lsun=sspLong, $ bsun=sspLat, $ NP=npole, $ Dg=distG, $ Dh=distH, $ typecoord='Apparent of date', $ timescale='UTC', $ resolution=img_reso, $ verfile=verfile, $ aperture=2.0, $ ;viewField=viewfield, $ omega0=2., $ epoch=epoch, $ objVertex=objVertex, $ objPolygons=objPolygons, $ albedoMap=albedoMap, $ shademod=shademod, $ ; onlyVR=1, $ image=image if (stars eq 1) && (s gt 1) && (counter eq 0) then begin backgrd = draw_stars(RA, DE, viewfield, image, LIST=list) endif image = image + backgrd if keyword_set(psf) then image = psf_convolve(image, psf=psf) image = bytscl(image, MAX=1.0, MIN=0.0) ;image = draw_label(image, img_reso, sepLong, sepLat, sspLong, sspLat, npole, targetName, epoch=epoch+0.0007*i, timescale) frames[*, *, i + counter] = image counter++ endif i = i + 1 endif endwhile FREE_LUN, unit ;------------------------------------- Saving the fits file and building the headers -------------------------------------; print, 'Saving file' print, outfile writefits, outfile, frames header = HEADFITS(outfile) ; SXADDPAR, header, 'BITPIX' , -32 SXADDPAR, header, 'OBJECT' , strtrim(targetName,2) ," Object's name" SXADDPAR, header, 'CRPIX1' , img_reso(0)/2. ,' (pix) Image X center' SXADDPAR, header, 'CRPIX2' , img_reso(1)/2. ,' (pix) Image Y center' SXADDPAR, header, 'CTYPE1' , 'RA---TAN' ,' Coordinate system of axis' SXADDPAR, header, 'CTYPE2' , 'DEC--TAN' ,' Coordinate system of axis' SXADDPAR, header, 'EQUINOX' , 2000. ,' Standard FK5' SXADDPAR, header, 'RADECSYS', 'FK5' ,' FK5' SXADDPAR, header, 'TYPECOOR', 'UTC' ,' Coordinates type' SXADDPAR, header, 'CD1_1' , -viewfield(0) / img_reso(0) ,' Translation matrix element' SXADDPAR, header, 'CD1_2' , 0 ,' Translation matrix element' SXADDPAR, header, 'CD2_1' , 0 ,' Translation matrix element' SXADDPAR, header, 'CD2_2' , viewfield(1) / img_reso(1) ,' Translation matrix element' SXADDPAR, header, 'CRVAL1' , crval1(0) ,' Right ascension of the fiducial point' SXADDPAR, header, 'CRVAL2' , crval2(0) ,' Declinaison of the fiducial point' ;SXADDPAR, header, 'DATE-OBS', ISOdate ,' Ephemeris date ('+timescale+')' ;SXADDPAR, header, 'SR#' , sr ,' Pole solution number (APDT)' ;SXADDPAR, header, 'L_SEP' , lsep(kDate) ,' (deg) Sub-Earth Point longitude' ;SXADDPAR, header, 'B_SEP' , bsep(kDate) ,' (deg) Sub-Earth Point latitude' ;SXADDPAR, header, 'L_SSP' , lsun(kDate) ,' (deg) Sub-Solar Point longitude' ;SXADDPAR, header, 'B_SSP' , bsun(kDate) ,' (deg) Sub-Solar Point latitude' ;SXADDPAR, header, 'NP' , NP(kDate) ,' (deg) Target North pole position angle' ;SXADDPAR, header, 'RA' , ra(kDate)/!DTOR ,' (deg) Object Right Ascension' ;SXADDPAR, header, 'DEC' , dec(kDate)/!DTOR ,' (deg) Object Declination' ;SXADDPAR, header, 'DIST_OBS', Dg(kDate) ,' (AU) Object range to observer' ;SXADDPAR, header, 'DSUN_OBS', Dh(kDate) ,' (AU) Object heliocentric distance' ;SXADDPAR, header, 'VMAG' , Vmag(kDate) ,' Object V magnitude' ;SXADDPAR, header, 'PHOTO_X' , photoCenter(0) ,' (pix) Photocentre X position' ;SXADDPAR, header, 'PHOTO_Y' , photoCenter(1) ,' (pix) Photocentre Y position' ;SXADDPAR, header, 'OMEGA' , omega0 ,' (rad/s) Object Spin' ;SXADDPAR, header, 'RMAX' , Rmax ,' (km) Object maximum radius' ;SXADDPAR, header, 'RAPP' , objScale(1)*2.0 ,' (mas) Object apparent radius' ;SXADDPAR, header, 'RADSHIFT', vmoy ,' (km/s) Observed mean radial velocity' ;SXADDPAR, header, 'VROT' , v2moy ,' (km/s) Rotational velocity integral (abs(vrad))' ;SXADDPAR, header, 'MUSOL' , musol(kDate) ,' (km/s) Sun radial velocity' ;SXADDPAR, header, 'MURHO' , murho(kDate) ,' (km/s) Earth radial velocity' ;SXADDPAR, header, 'MURA' , mura(kDate) ,' ("/s) Earth RA velocity' ;SXADDPAR, header, 'MUDEC' , mudec(kDate) ,' ("/s) Earth DEC velocity' ;SXADDPAR, header, 'KMPIX' , Rmax*winScale/240 ,' (km) Plate scale (km per pixel)' ;SXADDPAR, header, 'VERFILE' , locname ," Source file for Object's Shape (.ver)" ;SXADDPAR, header, 'OVERSAMP', oversamp ,' Pixels Oversampling' writefits, outfile, frames, header ; write_png, 'out.png', frames[*, *, 0] ;print, 'end' ;print, size(frames) END