;+ ; PRO VISU ; ; To obtain help, run: ; IDL> visu, /HELP ; ; To obtain informations on the development, run: ; IDL> visu, /INFO ; ;- pro visu, $ ;------------------------------------------------------------------------------------ ;-Object Orientation in Space lsep=lsep, $ ;-Sub-Earth Point Longitude (deg) bsep=bsep, $ ;-Sub-Earth Point Latitude (deg) lsun=lsun, $ ;-Sub-Solar Point Longitude (deg) bsun=bsun, $ ;-Sub-Solar Point Latitude (deg) NP=NP, $ ;-North to Pole Angle (NESW) (deg) ;-Object Ancilliary Information Vmag=Vmag, $ ;-Visual Apparent Magnitude omega0=omega0, $ ;-Target Sidereal Period Dg=Dg, Dh=Dh, $ ;-Geocentric and Heliocentric Distances (AU) RA=RA, DEC=DEC, $ ;-Target Equatorial Coordinates typecoord=typecoord, $ ;-Kind of Coordinates Information objectName=objectName, $ ;-String of the Object Name rMax=rMax, $ ;-Object Physical Size for Scale Display sr=sr, $ ;-Spin Solution as defined in APDT ;-Solar System Global State muSol=muSol, $ ;-Sun Radial Velocity muRho=muRho, $ ;-Earth Radial Velocity muRA=muRA, muDEC=muDEC, $ ;-Earth RA/DEC Velocity ;-Epoch of Observation, Timing Information epoch=epoch, $ timeScale=timeScale, $ ; Caracteristiques temporelles de l'objet ;------------------------------------------------------------------------------------ ;-Object Shape Model verFile=verFile, $ ;-Object Shape Model File ;-Object Display with Ellipsoid/Cellinoid Shape ellipsoid=ellipsoid, $ ;-Boolean to use Ellipsoid cellinoid=cellinoid, $ ;-Boolean to use Cellinoid ;-Object Display with Cellinoidal Shape axe_a1=axe_a1, axe_b1=axe_b1, axe_c1=axe_c1, $ ;-Axial Values for Ellipsoid & Cellinoid (Positive Side) axe_a2=axe_a2, axe_b2=axe_b2, axe_c2=axe_c2, $ ;-Axial Values for Cellinoid (Negative Side) ;-Diffusion Law DiffLawArr=DiffLawArr, $ ;-Relative Weight Lommel-Seeling vs Lambert ExpMinnaert=ExpMinnaert, $ ;-Minnaert Exponent shadeMod=shadeMod, $ ;------------------------------------------------------------------------------------ ;-Set Up Parameters for Display FITS=FITS, X=X, EPS=EPS, PNG=PNG, $ ;-Booleans to trigger X (default), FITS, EPS or PNG output ROTATION=ROTATION, $ ;-Boolean to Create an Animation of a Rotation outFile=outFile, $ ;-Filename upon Output (except X case) winScale=winScale, $ ;-Border between the Object and the Frame resolution=resolution, $ ;-X/Y-Axis Number of Pixel (default: [1024,768]) oversamp=oversamp, $ ;-Image Oversampling for FITS-Display Only labelINFO=labelINFO, $ viewField=viewField, $ ;------------------------------------------------------------------------------------ ;-Analysis Tools ;----------- ;-Keywords to be understood!!!!!!! ;----------- ;-Keywords to be developped fitellipse=fitellipse, $ ; VR=VR, onlyVR=onlyVR, $ ; Visualisation Option (cont'd) map_aster=map_aster, $ ; Visualisation ponderee par l'albedo maponly=maponly, $ ; Visualistion de l'albedo uniquement ;----------- ;-Keywords to be removed TBR mpeg=mpeg, $ ;--- >>> TBD change to ROTATION ;-???-; Rscal_x=Rscal_x, Rscal_y=Rscal_y, Rscal_z=Rscal_z, $ ; Object Characteristics, scaling of topographic model ; local=local, $ ;-RM RM RM;- a reflechir ; dirVER=dirVER, $ nom_aster=nom_aster, $ ;-TBR --> objectName pix=pix, $ ;-pixcel scle --- check to use a real fits with proper header image=image, $ ;------------------------------------------------------------------------------------ ;-Simulation of Observations from User's PSF PSF=PSF, $ ; To convolve the image with a theoritical PSF aperture=aperture, $ ; Telescope aperture (meter) wavelength=wavelength, $ ; Observed wavelength (micron) ;------------------------------------------------------------------------------------ ;-Developer's Corner INFO=INFO, $ ;-Display Author & History Information AUTHORS=AUTHORS, $ ;-Display Author Information HISTORY=HISTORY, $ ;-Display History Information PAUSE=PAUSE, $ ;-Twist Pause for X Outputs HELP=HELP, $ ;-Display Help/Howto on VISU objVertex=objVertex, $ objPolygons=objPolygons, $ nbVertex=nbVertex, $ nbFacet=nbFacet, $ albedoMap=albedoMap ;----------------------------------------------------------------------------------------------; ; ; ; 0 -- Developer's Corner, with Information and Help on the Software ; ; ; ;----------------------------------------------------------------------------------------------; ;-Print Everything if keyword_set( INFO ) then $ spawn, 'cat '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/header.info '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/history.info '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/authors.info ' ;-Print History of Modifications if keyword_set( HISTORY ) then $ spawn, 'cat '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/header.info '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/history.info ' ;-Print Author/Developer Names if keyword_set( AUTHORS ) then $ spawn, 'cat '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/header.info '+$ FILE_DIRNAME(ROUTINE_FILEPATH( 'visu' ))+'/docs/authors.info ' ;-Stop Program if keyword_set( INFO ) or $ keyword_set( HISTORY ) or $ keyword_set( AUTHORS ) then $ return ;----------------------------------------------------------------------------------------------; ; ; ; I -- Checking parameters, setting variables, reading files ; ; ; ;----------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------ ;-- I.1 -- Check & Set: Object Orientation in Space if not keyword_set(lsep) then lsep=0.D ;-Sub-Earth Point Longitude (deg) if not keyword_set(bsep) then bsep=0.D ;-Sub-Earth Point Latitude (deg) if not keyword_set(lsun) then lsun=0.D ;-Sub-Solar Point Longitude (deg) if not keyword_set(bsun) then bsun=0.D ;-Sub-Solar Point Latitude (deg) if not keyword_set(NP) then NP=0.D ;-North to Pole Angle (NESW) (deg) if not keyword_set(Dg ) then Dg=1.D ;-Geocentrique Distance (AU) if not keyword_set(Dh) then Dh=0.D ;-Heliocentrique Distance (AU) if not keyword_set(RA) then ra=0.D ;-Body Right Ascension (deg) if not keyword_set(DEC) then dec=0.D ;-Body Declination (deg) if not keyword_set(typeCoord) then $ ;-Type of Coordinates typeCoord='Apparent of date' ;------------------------------------------------------------------------------------ ;-- I.2 -- Check & Set: Object & Solar System Ancilliary Information if not keyword_set(muSol) then muSol=0.D ;-Sun Radial Velocity if not keyword_set(muRho) then muRho=0.D ;-Earth Radial Velocity if not keyword_set(muRA) then muRA =0.D ;-Earth RA Velocity if not keyword_set(muDEC) then muDEC=0.D ;-Earth DEC Velocity if not keyword_set(omega0) then omega0=0 ;-Object Rotation rate (rad/s) if not keyword_set(Vmag) then Vmag=0.D ;-Apparent Visual Magnitude if not keyword_set(sr) then sr=1 ;-Spin Solution Number from APDT ;------------------------------------------------------------------------------------ ;-- I.3 -- Check & Set: Epoch of Observation if keyword_set(epoch) then begin CALDAT, epoch, Month, Day, Year, Hour, Minute, Second ISOdate = string(Year, FORMAT='(%"%04I")')+'-'+$ string(Month, FORMAT='(%"%02I")')+'-'+$ string(Day, FORMAT='(%"%02I")')+'T'+$ string(Hour, FORMAT='(%"%2.2I")')+':'+$ string(Minute,FORMAT='(%"%2.2I")')+':'+$ string(Second,FORMAT='(%"%07.4F")') endif else begin ISOdate = '0000-00-00T00:00:00' endelse if not keyword_set(timeScale) then timeScale = 'UTC' ;------------------------------------------------------------------------------------ ;-- I.4 -- Check & Set: Display Parameters if not keyword_set(RESOLUTION) then $ ;-Image Dimensions (x,y) resolution=[1024,768] ;resolution=[1280,960] ;-Letter and Lines Look'n'Feel Definition letterSize = 1. * resolution(0)/640. letterThick = 1. * resolution(0)/640. lineThick = 1. * resolution(0)/640. ;-Default Device: Z-Buffer set_plot,'Z' ;-Set Up Display Device device, set_resolution=resolution printMod='' erase ;-Clean The Display Device if not keyword_set(pause) then pause=5 ;-Time Delay for X-Screen Display if not keyword_set(oversamp) then oversamp=1 ;-Pixels Oversampling (FITS mode only) if not keyword_set(winScale) then $ ;-Space between Object and Frame Borders winScale=1.2 ;-Labelling Information (Date, Orientation...) if not keyword_set(labelINFO) then begin labelINFO={name: 0, $ ;-Boolean to Print Object Name date: 0, $ ;-Boolean to Print Obsering Date orient:0, $ ;-Boolean to Print the Object Orientation (SEP, SSP, NP) frame: 0, $ ;-Boolean to Display the Celestial Reference Frame scale: 0, $ ;-Boolean to Display the Size Scale Bar text: [''], $ ;-User Labels (String Array) pos: [[0.9,0.1]] } ;-Position of User Label (NORMAL Coordinates) endif else begin ;-Look for already specified parameters tagNames=tag_names(labelINFO) nbTag=n_elements(tagNames) tagList=['NAME', 'DATE', 'FRAME', 'ORIENT', 'SCALE', 'TEXT', 'POS'] tagVals=[ 1, 1, 1, 1, 1, 1, 1 ] ;-Add non-specified tags to the structure for kTag=0, n_elements(tagList)-3 do begin tTag = where( tagNames eq tagList(kTag) ) if tTag(0) eq -1 then labelINFO=create_struct(tagList(kTag), tagVals(kTag), labelINFO) endfor tTag = where( tagNames eq 'TEXT' ) if tTag(0) eq -1 then labelINFO=create_struct('TEXT', [1.,1.], labelINFO) tTag = where( tagNames eq 'POS' ) if tTag(0) eq -1 then labelINFO=create_struct('POS', [1.,1.], labelINFO) endelse ;------------------------------------------------------------------------------------ ;-- I.5 -- Ellipsoids & Cellinoids -- Dimensions & Output Options dirVER = '/recherche/astrodata/shape_models/VER/' ;-TBR -- Shape models directory ;-Ellispoid (default); Nothing Specified if ( NOT keyword_set(verFile) and $ NOT keyword_set(cellinoid) and $ (NOT keyword_set(objVertex) OR $ NOT keyword_set(objPolygons))) then ellipsoid = 1 ;-Ellispoid (default): Model File & Cellinoid are Invalid if keyword_set(verFile) then $ if file_test(verFile) ne 1 then $ if NOT keyword_set(cellinoid) then ellipsoid = 1 ;-Creation of Ellipsoid or Cellinoid if keyword_set(ellipsoid) or keyword_set(cellinoid) then begin ;-Default Axes Dimensions if not keyword_set(axe_a1) then axe_a1= 1.2 if not keyword_set(axe_b1) then axe_b1= 1.0 if not keyword_set(axe_c1) then axe_c1= 0.8 ;-Axes Dimensions & I/O Names for Cellinoid if keyword_set(cellinoid) then begin if not keyword_set(axe_a2) then axe_a2 = 1.0 if not keyword_set(axe_b2) then axe_b2 = 0.8 if not keyword_set(axe_c2) then axe_c2 = 0.6 locName='cellinoid.ver' outfile='./cellinoid'+'' ;-Axes Dimensions & I/O Names for Ellipsoid endif else begin if not keyword_set(axe_a2) then axe_a2 = axe_a1 if not keyword_set(axe_b2) then axe_b2 = axe_b1 if not keyword_set(axe_c2) then axe_c2 = axe_c1 locName='ellipsoid.ver' outfile='./ellipsoid'+'' endelse ;-Create a Shape Model from Specified Dimensions create_ver_ellipsoid, objVertex, objPolygons, $ axes= [ [axe_a1, axe_b1, axe_c1], $ [axe_a2, axe_b2, axe_c2] ], $ nbVertex=nbVertex, $ nbFacet=nbFacet verFile = dirVER+locName endif ;-End of Test on Ellipsoid/Cellinoid ;------------------------------------------------------------------------------------ ;-- I.6 -- Load Shape Model ;-TBD ??? if keyword_set(verFile) then if file_test(verFile) eq 1 then begin locName = FILE_BASENAME(verFile, '.ver')+'.ver' ;-Read the Shape Model read_ver_shape_model, verFile, $ objVertex, objPolygons, $ nbVertex=nbVertex, $ nbFacet=nbFacet endif ;-Diffusion Law Defaults if not keyword_set(DiffLawArr) then DiffLawArr=[0.1,0.9,0.] ;-Relative Weigths Lommel-Seeliger vs Lambert if not keyword_set(ExpMinnaert) then ExpMinnaert=0.5 ;-Minnaert Coefficient for Minnaert's Diff. Law if not keyword_set(shadeMod) then shademod='polyshade' ;------------------------------------------------------------------------------------ ;-- I.6 -- Settings for Outputs if not keyword_set(outfile) then begin ;---- TBD - en attedant lien avec Eproc if keyword_set(objectName) then $ outfile='./'+objectName+'.'+printMod+'' $ else $ outfile='./'+'image_asteroide.'+printMod+'' endif ;-TBR-; if keyword_set(local) then outfile='./'+locName+'.'+printMod+'' else begin ;-TBR-; if keyword_set(nom_aster) then ima_asteroide=nom_aster else ima_asteroide='image_asteroide' ;---------------- keyword in the queue list if keyword_set(VR) then radialVelocity='VR' ;-radial velocity if keyword_set(onlyVR) then radialVelocity='onlyVR' ;-radial velocity if keyword_set(fitellipse) then mode_fitellipse='fitellipse' ;-ellipse analysis if keyword_set(map_aster) then mode_map='map_aster' ;-mapping if keyword_set(maponly) then mode_maponly='maponly' ;-mapping ;-Convolution PSF: diametre et lg onde if not keyword_set(aperture) then aperture = 0.2D0 if not keyword_set(wavelength) then wavelength = 0.5D0 ;-Object model scaling ;-???-; if not keyword_set(Rscal_x) then Rscal_x = 1.d0 ;-???-; if not keyword_set(Rscal_y) then Rscal_y = 1.d0 ;-???-; if not keyword_set(Rscal_z) then Rscal_z = 1.d0 ;---------------------------------------------------------------------------------------------------; ; Constant used in the soft ;--TBD -- quit it from hard dirMAP = '/astrodata/shape_models/VER/' ;-Albdeo maps directory ;-Alltime constants AUtoKM = 149597870.691d0 ;-Convertion AU/km DegToMAS = 3600000.d0 ;-Convertion degree/milli-arcsec textColor = 255 ; ;----------------------------------------------------------------------------------------------------; ; ;-- I.x -- Settings for Animations ; if keyword_set( ROTATION ) then begin ; ; ;-Set 360 Dates ; nbDate = 360 ; ; ;-Convering a Full 360 deg. Rotation ; lsep = lsep(0) + indgen(360) ; bsep = replicate( bsep(0), nbDate ) ; lsun = replicate( lsun(0), nbDate ) ; bsun = replicate( bsun(0), nbDate ) ; np = replicate( np(0), nbDate ) ; ; Dg = replicate( Dg(0), nbDate ) ; Dh = replicate( Dh(0), nbDate ) ; ; RA = replicate( RA(0), nbDate ) ; DEC = replicate( DEC(0), nbDate ) ; ; muRA = replicate( muRA(0), nbDate ) ; muDEC= replicate( muDEC(0), nbDate ) ; ; muSol = replicate( muSol(0), nbDate ) ; muRho = replicate( muRho(0), nbDate ) ; ; winScale = replicate( winScale(0), nbDate ) ; ; outfile='./tempo-'+string(indgen(nbDate), format='(I04)') ; ; endif ; kDate=0 ; ;-Loop over All the Date ; nbDate = n_elements( lsep ) ; for kDate=0, nbDate-1 do begin ;----------------------------------------------------------------------------------------------; ; ; ; II -- Shape model spatial orientation and facet properties ; ; ; ;----------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------ ;-- II.1 -- Current Frame Scale if not keyword_set(rMax) then rMax = max(sqrt( total(objVertex^2, 1) )) objScale = [ 0.5*rMax, (atan( .5*rMax / (Dg*AUtoKM) )/!DTOR)*DegToMAS ] ;print, objScale(1) ;-Facet Information Table FacetINFO = fltarr(15,nbFacet) if not keyword_set(albedoMap) then begin FacetINFO[13, *] = 1. endif else begin FacetINFO[13, *] = albedoMap endelse ;- 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) ;------------------------------------------------------------------------------------ ;-- II.2 -- Solar Direction & Spin Axis xSun = cos(lsun(kDate)*!DTOR)*cos(bsun(kDate)*!DTOR) ;-X ySun = sin(lsun(kDate)*!DTOR)*cos(bsun(kDate)*!DTOR) ;-Y zSun = sin(bsun(kDate)*!DTOR) ;-Z Dsun = [[xSun,ySun,zSun],[xSun,ysun,zSun]] ;-Distance SpinAxis = [[0d0,0d0,omega0],[0d0,0d0,omega0]] ;-In Object Reference Frame ;------------------------------------------------------------------------------------ ;-- II.3 -- Vertices Spatial Orientation T3D, /RESET ;- N T3D, rotate=[-90d0, -90d0, 0d0] ;-Permut Z -> Y -> X (IDL) E + W T3D, rotate=[0d0, -lsep(kDate), 0d0] ;-SEP Long at image center S T3D, rotate=[bsep(kDate), 0d0, 0d0] ;-SEP Lat at image center T3D, rotate=[0d0, 0d0, NP(kDate)] ;-North to Pole Angle alignment ;-- Recuperation of pole axis vector paxis=[0., 0., 1., 1.] paxis = !P.T ## paxis ;------------------------------------------------------------------------------------ ;-- II.4 -- Rotate the shape model obsVertex = Vert_T3D( objVertex ) SpinAxis = Vert_T3D( SpinAxis, double=1 ) Dsun = Vert_T3D( Dsun, double=1 ) Dsun = Dsun(*,0) Dearth = [0,0,1] ;-SEP is along Z in IDL reference frame SpinAxis = SpinAxis[*,0] ;---TBD if keyword_set(viewField) then begin R = [Dg*AUtoKM*tan(viewField(0)*!DTOR/2), Dg*AUtoKM*tan(viewField(1)*!DTOR/2)] endif else begin R = [winScale(kDate)*rMax, winScale(kDate)*rMax] endelse ;------------------------------------------------------------------------------------ ;-- III.1 -- Scale and Shading Direction (Light from the Sun...) scale3, xrange=float(resolution(0))/resolution(1)*[-R(0),R(0)], $ yrange=[-R(1),R(1)], zrange=[-max(R), max(R)] ;---------------------------- Polyshade -------------------------------------------- ;lum=long(1000/((Dg+Dh)^2)) lum = 255 set_shading, values=[2,lum], LIGHT=Dsun, $ ;-Coordinate West-North-Earth /reject ;-- Draw the model shaded with values of albedo mapShade = fltarr(2048, 2048) + 1 mapShade = POLYSHADE(obsVertex, objPolygons, $ POLY_SHADES=FacetINFO[13, *], $ /T3D) mapShade = tvrd() mapShade = float(mapShade / 255.) ERASE if shadeMod eq 'polyshade' then begin ;-- Draw the model shaded as lighten from the sun lightShade = fltarr(2048, 2048) + 1 lightShade = POLYSHADE(obsVertex, objPolygons, $ /T3D) lightShade = tvrd() lightShade = float(lightShade / 255.) ;-- combining both shadings (light * albedo) image = mapShade * lightShade ;plot, float(resolution(0))/resolution(1)*[-R,R], [-R,R], /NORMAL, /NOERASE, $ ; position= [0,0,1,1], $ ; xst=5, yst=5, $ ; /NODATA ; ;-- Poles axis ; XX = 2*Rmax*[-paxis(0), paxis(0)] ; YY = 2*Rmax*[-paxis(1), paxis(1)] ; oplot, XX, YY, linestyle=0, thick=lineThick ;----TBR ; ;-- III.4.2 -- Object Orientation on the Plane of the Sky ; if keyword_set( labelINFO.orient ) then begin ; ; ;-Sub-Earth Point Coordinates ; XYouts, 0.035, 0.105, /NORMAL, $ ; 'SEP !4k!3,!4b!3 = '+string (lsep(kDate), FORMAT='(%"%5.1f")')+'!Uo!N'+$ ; ','+string (bsep(kDate), FORMAT='(%"%5.1f")')+'!Uo!N', $ ; charsize=1.2*letterSize, $ ; charthick=1.2*letterThick ; ; ;-Sub-Solar Point Coordinates ; XYouts, 0.035, 0.065, /NORMAL, $ ; 'SSP !4k!3,!4b!3 = '+string (lsun(kDate), FORMAT='(%"%5.1f")')+'!Uo!N'+$ ; ','+string (bsun(kDate), FORMAT='(%"%5.1f")')+'!Uo!N', $ ; charsize=1.2*letterSize, $ ; charthick=1.0*letterThick ; ; ;-Pole to North Angle ; XYouts, 0.035, 0.025, /NORMAL, $ ; 'NP = '+string (NP(kDate), FORMAT='(%"%5.1f")')+'!Uo!N', $ ; charsize=1.2*letterSize, $ ; charthick=1.0*letterThick ; endif ; ; ; ;-- III.4.3 -- Body's Name in Uppercase Letters ; if keyword_set( labelINFO.name ) then $ ; XYouts, 0.035, 0.15, /NORMAL, $ ; charsize=1.75*letterSize, $ ; charthick=1.5*letterThick, $ ; '!5'+strupcase(FILE_BASENAME(locName,'.ver')) ; ; ;-- III.4.4 -- Observing Date ; if keyword_set( labelINFO.date ) then $ ; XYouts, 0.035, 0.95, /NORMAL, $ ; strmid(string(ISODATE),0,21)+ ' '+timescale, $ ; charsize=1.5*letterSize, $ ; charthick=1.5*letterThick ; ; ; ;-- III.4.5 -- Celestial Frame ; if keyword_set( labelINFO.name ) then begin ; ; ;-Horizontal Axis ; XX = [-R, R] ; YY = [ 0, 0] ; oplot, XX, YY, linestyle=2, thick=lineThick ; ; ;-Vertical Axis ; XX = [ 0, 0] ; YY = [-R, R] / 1.2 ; oplot, XX, YY, linestyle=2, thick=lineThick ; ; ;-Celestial North and East ; XYouts, -R*1.05, -0.02*R, 'E', charsize=1.2*letterSize, charthick=1.25*letterThick, /DATA ; XYouts, 0.0, R/1.2, 'N', charsize=1.2*letterSize, charthick=1.25*letterThick, /DATA, align=0.5 ; ; endif ; ; ;-- III.4.6 -- Scale Segment Caption (km and apparent size) ; if keyword_set( labelINFO.scale ) then begin ; ; ;-Drawing the Scale Segment ; YY = [-R*.94,-R*.94] ; XX = [0,(Rmax/2)] ; oplot, XX, YY, thick=lineThick ; ; ;-Print size in KM and MAS ; xyouts, 0, -.92*R, /DATA, $ ; strtrim(string(objScale(0), FORMAT='(F9.2)'),2)+' km '+'- '+$ ; strtrim(string(objScale(1), FORMAT='(F9.2)'),2)+' mas', $ ; charsize=.7*letterSize, charthick=.6*letterThick ; ; endif ; ; ; ; XX = [-R, R] ; YY = [ 0, 0] ; oplot, XX, YY, linestyle=2, thick=lineThick ; ;----TBR ; image = tvrd() endif else if shadeMod eq 'lsl' then begin ;-IV.1. Vertices Cartesian Coordinate (for each Facet) vertX = fltarr(3, NbFacet) ;-X Coordinate vertY = fltarr(3, NbFacet) ;-Y Coordinate vertZ = fltarr(3, NbFacet) ;-Z Coordinate for kFacet=0L, NbFacet-1L do begin FacetINFO[0:2, kFacet] = objPolygons[ 4*kFacet+1:4*kFacet+3 ] vertX[*, kFacet] = objVertex[ 0, objPolygons[4*kFacet+1:4*kFacet+3] ] vertY[*, kFacet] = objVertex[ 1, objPolygons[4*kFacet+1:4*kFacet+3] ] vertZ[*, kFacet] = objVertex[ 2, objPolygons[4*kFacet+1:4*kFacet+3] ] endfor ;-IV.2. Facet Surface's Normal Vector Mesh_normals = compute_mesh_normals2( obsVertex(0:2,*), objPolygons ) DotProduct_Sun = Mesh_normals##transpose(Dsun) DotProduct_Earth = Mesh_normals##transpose(Dearth) ;-IV.3. Illuminated Faces ;visibleFacet = where( (DotProduct_Sun gt 0 and DotProduct_Earth gt 0), nbIllum ) visibleFacet = where( (DotProduct_Earth gt 0), nbIllum ) FacetINFO( 3, visibleFacet ) = 1 illumVertex = intarr(3*nbIllum) for kIllum=0, nbIllum-1 do illumVertex(3*kIllum:3*kIllum+2) = FacetINFO(0:2, visibleFacet(kIllum)) illumVertex = illumVertex( uniq( illumVertex, sort(illumVertex) ) ) ;-IV.4. Facets with Shadow (Concavities or Binary Bodies) T3D, /RESET T3D, ROTATE = [-90d0, -90d0, 0d0] T3D, ROTATE = [0d0, -lsun(kDate), 0d0] T3D, ROTATE = [bsun(kDate), 0d0, 0d0] ;-Shape Model as seen from the Sun verTest = Vert_T3D(objVertex) ;-Compute Barycenter of each Facet baryArr=fltarr(3,nbIllum) for kFacet=0, nbIllum-1 do $ baryArr[*, kFacet] = total( verTest[*, objPolygons[4 * visibleFacet[kFacet] + [1, 2, 3]]], 2 ) / 3. ;-Look for Facets one in front of another for kFacet=0, nbIllum-1 do begin ;-Distance between the Barycenters distance = sqrt( (transpose(baryArr[0, *] - baryArr[0, kFacet]))^2. + $ (transpose(baryArr[1, *] - baryArr[1, kFacet]))^2. ) ;-Sorting the Barycenters per Distance Order ordreDIST=sort( distance ) ;-Finding out How Many Facets should be tested (min=optimized code) bordVert = verTest[*, objPolygons[4 * visibleFacet[kFacet] + [1, 2, 3]]] bordVert[0, *] = ( bordVert[0, *] - baryArr[0, kFacet] )^2. bordVert[1, *] = ( bordVert[1, *] - baryArr[1, kFacet] )^2. bordVert[2, *] = ( bordVert[2, *] - baryArr[2, kFacet] )^2. bordDist = sqrt( total( bordVert, 1) ) trash=where( distance le bordDist, nbToCheck) ;-Looking for any aligned Facets (seen from the Sun) for kBary=1, nbToCheck-1 do begin ;-Looking for any aligned Facets (seen from the Sun) ;-Retrieve the Cartesian Coordinates of the Vertexes Associated with the Barycenter cornerArr = objPolygons[4 * visibleFacet[ordreDIST[kBary]] + [1, 2, 3]] X1 = verTest[0, cornerArr[0]] & X2 = verTest[0, cornerArr[1]] & X3 = verTest[0, cornerArr[2]] Y1 = verTest[1, cornerArr[0]] & Y2 = verTest[1, cornerArr[1]] & Y3 = verTest[1, cornerArr[2]] Z1 = verTest[2, cornerArr[0]] & Z2 = verTest[2, cornerArr[1]] & Z3 = verTest[2, cornerArr[2]] ;-Facet triangle delimiting lines slope12= ( Y2-Y1 ) / ( X2-X1 ) slope13= ( Y3-Y1 ) / ( X3-X1 ) slope23= ( Y3-Y2 ) / ( X3-X2 ) ordOrigin12= Y1 - slope12*X1 ordOrigin13= Y1 - slope13*X1 ordOrigin23= Y2 - slope23*X2 ;-Checking if the Barycenter is inside the Facet if ( (baryArr[1, kFacet] - slope12 * baryArr[0, kFacet] - OrdOrigin12) * (Y3 - slope12 * X3 - OrdOrigin12) gt 0 ) and $ ( (baryArr[1, kFacet] - slope13 * baryArr[0, kFacet] - OrdOrigin13) * (Y2 - slope13 * X2 - OrdOrigin13) gt 0 ) and $ ( (baryArr[1, kFacet] - slope23 * baryArr[0, kFacet] - OrdOrigin23) * (Y1 - slope23 * X1 - OrdOrigin23) gt 0 ) $ then begin ;-Selecting the Facet BEHIND to be shadowed if baryArr[2, kFacet] ge baryArr[2, ordreDIST[kBary]] then $ facetINFO[3, visibleFacet[ordreDIST[kBary]]] = 0 else $ facetINFO[3, visibleFacet[kFacet]] = 0 ;-Once the Facet already shadowed, interrupt the loop goto, nextBarycenter endif endfor ;-End of loop over barycenters nextBarycenter: endfor ;-End of loop over Facets ;-Retrieve the Illuminated Facets, Taking the Shadows into Account visibleFacet = where( FacetINFO[3, *] eq 1, nbIllum ) shadowedFacet = where( FacetINFO[3, *] eq 0, nbShadowed ) ;-Sorting Facet from the Background to the Observer zMax = max( vertZ, dimension=1) zOrder = sort( zMax ) zOrderVis = sort( zMax[visibleFacet] ) zOrderSha = sort( zMax[shadowedFacet]) ;-IV.4. Reflection (4) and Incidence (5) Angles FacetINFO(4,*) = DotProduct_Earth FacetINFO(5,*) = DotProduct_Sun ;-IV.5. Diffusion Laws FacetINFO( 9,*)= FacetINFO(5,*) ;-Lambert's Law FacetINFO(10,*)= FacetINFO(5,*)/( FacetINFO(4,*)+FacetINFO(5,*) ) ;-Lommel-Seelinger's Law FacetINFO(11,*)= ( FacetINFO(5,*)^(ExpMinnaert) )*( FacetINFO(4,*)^(ExpMinnaert-1) ) ;-Minnaert's Law ;-IV.6. Radial Velocity Computation for kFacet=0L, nbFacet-1L do begin ;-Radial Velocity Computation vecFacet = [ mean(vertX(*,kFacet)), mean(vertY(*,(kFacet))), mean(vertZ(*,kFacet)) ] FacetINFO(6,kFacet)=-TRANSPOSE(Dearth)#CROSSP(SpinAxis, vecFacet ) FacetINFO(8,kFacet)=-TRANSPOSE(Dsun ) #CROSSP(SpinAxis, vecFacet ) ;-Facet Surface FacetINFO(7, kFacet)= Poly_area( vertX(*, kFacet), vertY(*, kFacet) ) endfor ;-IV.7. Intensity for each Facet FacetINFO(12, *) = FacetINFO(7, *) * FacetINFO(4, *) * ( DifflawArr(0) * FacetINFO( 9, *) $ + DifflawArr(1) * FacetINFO(10, *) $ + DiffLawArr(2) * FacetINFO(11, *) ) ;-Read again the model, this time as seen from the Earth vertX = fltarr(3,nbFacet) ;-X Coordinate vertY = fltarr(3,nbFacet) ;-Y Coordinate vertZ = fltarr(3,nbFacet) ;-Z Coordinate for kFacet=0L, nbFacet-1L do begin vertX(*,kFacet) = obsVertex( 0, objPolygons(4*kFacet+1:4*kFacet+3) ) vertY(*,kFacet) = obsVertex( 1, objPolygons(4*kFacet+1:4*kFacet+3) ) vertZ(*,kFacet) = obsVertex( 2, objPolygons(4*kFacet+1:4*kFacet+3) ) endfor ;-IV.8. Total Surface and Total Intensity intensityArr=fltarr(2,2) intensityArr(0,0) = total( FacetINFO(7,visibleFacet) ) ;-Total Visible Surface intensityArr(1,0) = total( FacetINFO(12,visibleFacet) ) ;-Total Visible Intensity intensityArr(0,1) = total( FacetINFO(7,*) ) ;-Total Surface intensityArr(1,1) = total( FacetINFO(12,visibleFacet) ) ;-Total Intensity ; polyfill, vertX[*, shadowedFacet[zOrderSha[kFacet]]], $ ;---- TBD -- A quoi correspond quoi??? ;-IV.9. Speed & Square-Speed Computation speedArr = fltarr(3,2) speedArr(0,0) = total( ( FacetINFO( 6,visibleFacet) + FacetINFO( 8,visibleFacet) ) $ ;-Mean Visible Speed * FacetINFO(12,visibleFacet) * FacetINFO(13,visibleFacet) ) / intensityArr(1,0) speedArr(1,0) = total( abs( (FacetINFO( 6,visibleFacet) + FacetINFO( 8,visibleFacet)) $ * FacetINFO(12,visibleFacet) * FacetINFO(13,visibleFacet) ) ) / intensityArr(1,0) ;-Mean Visible Square Speed speedArr(2,0) = max( abs(FacetINFO(6,visibleFacet)+FacetINFO(8,visibleFacet)) ) ;-Max Visible Speed speedArr(0,1) = total( ( FacetINFO( 6,*) + FacetINFO( 8,*) ) $ ;-Mean Speed * FacetINFO(12,*) * FacetINFO(13,*) ) / intensityArr(1,1) speedArr(1,1) = total( abs( (FacetINFO( 6,*) + FacetINFO( 8,*)) $ * FacetINFO(12,*) * FacetINFO(13,*) ) ) / intensityArr(1,1) ;-Mean Square Speed speedArr(2,1) = max( abs(FacetINFO(6,*)+FacetINFO(8,*)) ) ;-Max Speed ;-IV.10. Printing the Radial Velocity Model if keyword_set(radialVelocity) then begin; print, 'radial velocity' ;-Erase the Z-buffer polyfill,[0,0,1,1],[0,1,1,0], /IMAGE_INTERP,PAT=[127,127,127] ;-Radial Velocity Frame Creation for kFacet=0L, nbFacet-1L do $ polyfill, .5+vertX[*,zOrder(kFacet)]/(2*R)*480/640, $ .5+vertY[*,zOrder(kFacet)]/(2*R), /IMAGE_INTERP, $ PAT=[127,127,127]+( FacetINFO(6,zOrder(kFacet))+FacetINFO(8,zOrder(kFacet)))/speedArr(2)*127. ;-Write the Frame in a variable if printMod eq 'fits' then imageVR=tvrd() ;-Erase the Z-buffer polyfill,[0,0,1,1],[0,1,1,0],/IMAGE_INTERP,PAT=[0,0,0] endif ;---TBR vmoy= speedArr(0,0) v2moy= speedArr(1,0) ;---- tout ca est inutile. Soit garder speedArr, soit juste creer vmo v2moy ; vmoy= speedArr(0,1) ; v2moy= speedArr(1,1) ;---TBR ;print, vmoy, v2moy ;print, speedArr FacetINFO[13, *] = FacetINFO[13, *] / 255. ; -IV.11. Printing the Model with Diffusion Law ; for kFacet=0L, nbShadowed-1L do begin ; polyfill, vertX[*, shadowedFacet[zOrderSha[kFacet]]], $ ; vertY[*, shadowedFacet[zOrderSha[kFacet]]], $ ; /IMAGE_INTERP, $ ; PAT = [255, 255, 255] ; endfor for kFacet=0L, nbIllum-1L do begin polyfill, vertX[*, visibleFacet[zOrderVis[kFacet]]], $ vertY[*, visibleFacet[zOrderVis[kFacet]]], $ /IMAGE_INTERP, $ PATTERN = [255, 255, 255] * FacetINFO[ 4, visibleFacet[zOrderVis[kFacet]]] $ * FacetINFO[13, visibleFacet[zOrderVis[kFacet]]] *$ ( DiffLawArr[0] * FacetINFO[10, visibleFacet[zOrderVis[kFacet]]] $ + DiffLawArr[1] * FacetINFO[11, visibleFacet[zOrderVis[kFacet]]] $ + DiffLawArr[2] * FacetINFO[12, visibleFacet[zOrderVis[kFacet]]] ) endfor ; endelse ; 6. Calcul de l'etallement de la vitesse : variance ; v = 0 ; sigmav = 0 ; vMax= speedArr(2,0) ; surface_visible = intensityArr(0,0) ;-to be removed ; Itotal = intensityArr(1,0) ;-to be removed ; Itotal =intensityArr(0,1) ; vMax_2 =speedArr(2,1) ; surface_visible_2=intensityArr(1,1) ; ;gv=fltarr(2,20) ; for k=0L,nbFacet-1L do begin ; l = zOrder[k] ; if FacetINFO[3,FacetINFO[6,l]] eq 1 then begin ; sigmav=sigmav+(vmoy-FacetINFO[7,l])^2*FacetINFO[8,l]*(FacetINFO[4,l]*(DiffLawArr(1)*FacetINFO[11,l])+DiffLawArr(0)*FacetINFO[5,l])/Itotal ; endif ; endfor lightShade = tvrd(/ORDER) lightShade = float(lightShade / 255.) ;image = mapShade * lightShade image = lightShade endif ; ;------------------------------------------------------------------------------------ ; ;-- V.3 -- Compute the Position of Photocenter ; photoCenter = [ total( total(image,2) * indgen(resolution(0)) ), $ ; total( total(image,1) * indgen(resolution(1)) ) ] / total(image) ; ;;--- TBD---- ; if keyword_set(mode_fitellipse) then begin ; mask = where( image ne 0) ; imSize = size(image) ;; ellipsePts = fit_ellipse(mask, xsize=imSize(1), ysize=imSize(2), axes=axes, orientation=orientation) ; endif ;------------------------------------------------------------------------------------ ; ;-- V.4 -- Radial Velocities ; if keyword_set(radialVelocity) then begin ; ; ;-V.4.1. Cube Creation: 1-Photometry, Radial Velocity, Radial Velocity^2 ; if radialVelocity eq 'VR' then begin ; ;-Cube Creation ; cube = fltarr( resolution(0), resolution(1), 3) ; cube(*,*,0) = image ;-Image ; cube(*,*,2) = speedArr(2,1)*(imageVR-127)/127 ;-VR, all Facets ; ; ;-VR, visible Facets ONLY ; maskBody = where(image eq 0) ; imageVR2 = speedArr(2,0)*(imageVR-127)/127 ; imageVR2(maskBody) = 0 ; cube(*,*,1) = imageVR2 ; ; image=cube ; endif ; ; ;-V.4.2. Radial Velocity ONLY ; if radialVelocity eq 'onlyVR' then image = speedArr(2,0)*(imageVR-127)/127 ; ; endif ;-V.5. Image OverSampling ; if keyword_set(pix) then begin ; ; ;-pix = Nb pixel for 2*Rmax ; pix = oversamp*pix ; image = tvrd() ; ; scaleY = indgen(floor(winScale(kDate)*pix)) *480L/(winScale(kDate)*pix) ; scaleX = indgen(floor(winScale(kDate)*pix*640/480))*480L/(winScale(kDate)*pix) ; ; image = interpolate(image,scaleX,scaleY,/GRID) ; ; endif ; ; ;-V.6. Write the image WITH HEADER ; writefits, outfile, image, header ; ; ; ; endif else begin ;-If NOT a FITS file ; ; ;-V.6. Overplotting an Ellipse ; if keyword_set(mode_fitellipse) then begin ; ; ;-Recover image ; image = tvrd() ; mask = where( image ne 0) ; ; ;-Body Apparent Ellipse ; imSize = size(image) ;; ellipsePts = fit_ellipse(mask, xsize=imSize(1), ysize=imSize(2), axes=axes, orientation=orientation) ; plots,ellipsepts(*,*),/device, Color=FSC_Color('red') ; XYouts,400,50,charsize=1,'2a,2b = '+string(axes(0)*objScale(1)/(120/winScale(kDate)),FORMAT='(%"%7.1f")')+', ' $ ; +string(axes(1)*objScale(1)/(120/winScale(kDate)),FORMAT='(%"%7.1f")')+' arcsec', /DEVICE ; ; endif ; ;-Recover image ; image = tvrd() ; ; ;-Body's Name in Uppercase Letters ; XYouts, 30, 80,charsize=2,'!5'+strupcase(FILE_BASENAME(locname,'.ver')),/DEVICE ; ; ;-Sub-Earth Point & Sub-Solar Point Coordinates and Pole Angle ; XYouts, 30, 55,'!4k!3 SEP='+string(lsep(kDate), FORMAT='(%"%5.1f")')+'!3�'+$ ; ' !4b!3 SEP='+string(bsep(kDate), FORMAT='(%"%5.1f")')+'!3�', /DEVICE ;-SEP ; XYouts, 30, 35,'!4k!3 SSP='+string(lsun(kDate), FORMAT='(%"%5.1f")')+'!3�'+$ ; ' !4b!3 SSP='+string(bsun(kDate), FORMAT='(%"%5.1f")')+'!3�', /DEVICE ;-SSP ; XYouts, 30, 15,'NP='+string(NP(kDate), FORMAT='(%"%5.1f")')+'!3�', /DEVICE ;-Pole Angle ; ; ;-Celestial Frame ; XYouts, 30, 250,charsize=2,'E',charthick=2, /DEVICE ;-East ; XYouts, 315, 470,charsize=2,'N',charthick=2, /DEVICE ;-North ; ;-Horizontal Axis ; YY=[.5,.5] ; XX=[.1,.9] ; oplot,XX,YY,linestyle=2 ; ;-Vertical Axis ; YY=[.1,.9] ; XX=[.5,.5] ; oplot,XX,YY,linestyle=2 ; ; ;-Drawing the Scale Segment ; YY=[.03,.03] ; XX=[.5,(320+(240/winScale(kDate))/2)/640] ; oplot,XX,YY,thick=2 ; Affichage de la barre d'echelle ; ; ;-Scale Segment Caption (km and apparent size) ; XYouts,640/2,23,string(objScale(0), FORMAT='(%"%5.1f")')+'km'+' - '+string(objScale(1), FORMAT='(%"%9.1f")')+'mas', /DEVICE ; endelse ; endfor ;-End of loop over dates end