;+ ; NAME: ; RemoteView_FOV_Cube ; PURPOSE: ; CATEGORY: ; 3D viewing ; CALLING SEQUENCE: FUNCTION RemoteView_FOV_Cube, FOVInfo, x, R=R, bcube=BCube ; INPUTS: ; FOVInfo array[1]; type: structure ; OUTPUTS: ; status 0: no data inside field of view for specified viewing direction ; 1: there are data inside the field of view ; FOVCube array[nX,nY,nZ]; type: float ; function values covering fov of nX,nY lines of sight ; at nZ positions along line of sight ; OPTIONAL OUTPUT PARAMETERS: ; R=R array[3,nX,nY,nZ]; type: float ; heliographic locations at all of the line of sight segments. ; Always returned in spherical coordinates (even for rectangular ; matrix) i.e. heliographic longitude, latitude in radians; and ; heliocentric distance in AU. ; bcube=BCube array[nX,nY,nZ]; type: float ; function values for magnetic field covering fov of nX,nY lines of sight ; at nZ positions along line of sight ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; sphere_distance, cvt3d, RemoteView_FOV_xyz, AngleRange ; PROCEDURE: ; > Read the RemoteView document (in D:\work\doc\remoteview.doc) for a detailed description. ; > The 3D data set is specified in an arbitrary Cartesian coordinate ; system x-y-z. ; > The data set is contained within a sphere of radius DataRadius centered ; on the origin. ; > Put a rectangle with sides 2*tan(FOVInfo.view_fov) in a plane perpendicular ; to the viewer direction at a distance of one unit. Define a rectangular grid ; of dimension nX x nY covering the rectangle. Connect the origin with each of ; the nX x nY grid points to define an angular grid. ; > Define a radial grid of dimension nZ extending outward from the ; observer between FOVInfo.view_extent[0] and FOVInfo.view_extent[1]. ; > The radial and angular grid define a 'sky grid' covering the portion of ; space along the viewing direction which may contain data. ; MODIFICATION HISTORY: ; FEB-1998, Paul Hick (UCSD/CASS) ; DEC-1998, Paul Hick (UCSD/CASS, pphick@ucsd.edu), rewrite ;- rInner = FOVInfo.matrix_rEdge[0] rOuter = FOVInfo.matrix_rEdge[1] Loc = FOVInfo.view_loc Dir = FOVInfo.view_dir ; The calculation of the sky grid starts here. First check whether for the specified ; viewing direction the data sphere lies inside the field of view. If not, then ; return with Status=0 IF Loc[2] GE rOuter THEN BEGIN ; If viewer outside data sphere ; The whole square field of view of size 2*Fov is included in a cone with ; opening angle atan(sqrt(2.)*tan(Fov). For a viewer outside the data ; sphere the data fit in a cone with opening angle asin(rOuter/Loc[2]). ; If the elongation of the viewing direction is greater than the sum of the two ; cone angles then the fov does not intersect the data sphere. ; tan(fov) is the linear half-width of the fov in plane perpendicular to the ; viewing direction at unit distance from viewing location fovlim = tan(FOVInfo.view_fov) fovlim = atan( sqrt( total(fovlim*fovlim) ) ) + asin(rOuter/Loc[2]) ;fovlim = atan(sqrt(2.)*tan(max(FOVInfo.view_fov)))+asin(rOuter/Loc[2]) IF !pi-sphere_distance(Loc[0:1], Dir) GT fovlim THEN BEGIN message, /info, 'data sphere outside field of view' RETURN, 0 ENDIF ENDIF nDim = FOVInfo.view_ndim nDimX = nDim[0] nDimY = nDim[1] nDimZ = nDim[2] x = fltarr(4,nDimX,nDimY,nDimZ, /nozero); Don't change this! x = RemoteView_FOV_xyz(x, FOVInfo) ; Fill x[0:2,*,*,*] with x,y,z coordinates x[3,*] = 1. ; Needed for transformation with pt ; x now contains rectangular coordinates in the X-Y-Z coordinate system (with viewer in the ; origin, Z along the view direction and X and Y along the sides of the field of view). ; First we need the elongations across the field of view. ; Calculate location of Sun in the X-Y-Z coordinate system ;sLoc = [!pi+Loc[0],-Loc[1],Loc[2]] ;sLoc = EulerRotate([!pi+Dir[0],-(!pi/2.-Dir[1]),-FOVInfo.view_tilt], sLoc) ;sLoc = cv_coord(from_sphere=sLoc, /to_rect) ; Calculate elongation using the depth layer in x which is farthest away from the viewer location ;*FOVInfo.view_elo = sphere_distance(x[0:2,0:nDimX*nDimY-1], sLoc) ; Transform the fov coordinates to the original x-y-z coordinates. ; Note: a rotation of a coordinate system means applying the opposite rotation ; to the vector-components. Same for translation. ; 1. Rotate fov coordinates from X-Y,Z to x"-y"-z" (rotation over -Tilt) ; 2. Rotate x"-y"-z" to x'-y'-z'. The Euler angles for this rotation are ; Alfa=0., Beta=!pi/2.-Dir[1], Gamma=!pi-Dir[0] ; 3. Translate over -Loc from x'-y'-z' to x-y-z ; 4. Rotate over the heliographic longitude offset lng0 Loc = cv_coord(from_sphere=Loc, /to_rect) mat = cvt3d(rot=[0.,0.,-FOVInfo.view_tilt, 0.,!pi/2.-Dir[1],!pi-Dir[0]]) mat = cvt3d(mat=mat, trans=-Loc, rot=[0.,0.,FOVInfo.matrix_lng0]) ; The matrix 'mat' will convert x to heliocentric, rectangular coordinates. ; The x-axis is lined up with the matrix start. The distance unit is AU. R = mat#x ; Rectangular heliocentric coordinates (AU) ;============================================================================= ; Do-loop is faster in low memory configuration ; ;R = reform(R,4,nDimX*nDimY,nDimZ, /overwrite) ;for i=0,nDimZ-1 do R[0:2,*,i] = cv_coord(from_rect=R[0:2,*,i], /to_sphere) ;R = reform(R,4,nDimX*nDimY*nDimZ, /overwrite) ; ;============================================================================= R[0:2,*] = cv_coord(from_rect=R[0:2,*], /to_sphere) ; Spherical heliocentric coordinates (AU) R[0,*] = AngleRange(R[0,*]) ; Longitude must be in [0,2*!pi] CASE FOVInfo.matrix_rect OF ; If interpolating on rectangular data 0: BEGIN ; Scale to units of array index. The translation shifts the latitude values from ; [-90,+90] to [0,180], and shifts the heliocentric distance by rInner (AU). ; The scaling converts from degrees and AU to grid spacings. sz = size(*FOVInfo.matrix, /dim) ; Three spatial dimensions mat = cvt3d(trans=[0.,-!pi/2,rInner], scale=[2*!pi/(sz[0]-1),!pi/(sz[1]-1),FOVInfo.matrix_dR]) x = mat#R END 1: BEGIN ; Now convert to array indices which can be used to interpolate on the the data array. mat = cvt3d(mat=mat, scale=FOVInfo.matrix_dR) ; Convert distance scale to grid spacings mat = cvt3d(mat=mat, trans=-FOVInfo.matrix_origin) ; Translate over -origin. x = mat#x ; Convert to array indices END ENDCASE ;print, 'new' ;print, min(x[0,*]),max(x[0,*]) ;print, min(x[1,*]),max(x[1,*]) ;print, min(x[2,*]),max(x[2,*]) IF FOVInfo.matrixb_index NE -1 THEN BEGIN BCube = interpolate((*FOVInfo.matrixb)[*,*,*,FOVInfo.matrixb_index], x[0,*],x[1,*],x[2,*], missing=0.) BCube = reform(BCube, nDim, /overwrite) ENDIF x = interpolate((*FOVInfo.matrix)[*,*,*,FOVInfo.matrix_index], x[0,*],x[1,*],x[2,*], missing=0.) x = reform(x, nDim, /overwrite) ; Ouput nDimX x nDimY x nDimZ array R[0,*] = AngleRange(FOVInfo.matrix_lng0+R[0,*]) ; Back to heliographic longitudes R = reform(R[0:2,*], [3,nDim]) RETURN, 1 & END