[cig-commits] r6656 - geodyn/3D/MAG/trunk/idl
wei at geodynamics.org
wei at geodynamics.org
Tue Apr 24 12:55:01 PDT 2007
Author: wei
Date: 2007-04-24 12:55:01 -0700 (Tue, 24 Apr 2007)
New Revision: 6656
Added:
geodyn/3D/MAG/trunk/idl/magmovCMB.pro
geodyn/3D/MAG/trunk/idl/magmovEQ3.pro
Log:
Add two more postprocessing programs, magmovCMB displays composite images of Br on the outer spherical surface in three map view: aitoff, north and south polar. magmovEQ3 displays results in the equatorial plane.
Added: geodyn/3D/MAG/trunk/idl/magmovCMB.pro
===================================================================
--- geodyn/3D/MAG/trunk/idl/magmovCMB.pro 2007-04-24 19:40:51 UTC (rev 6655)
+++ geodyn/3D/MAG/trunk/idl/magmovCMB.pro 2007-04-24 19:55:01 UTC (rev 6656)
@@ -0,0 +1,439 @@
+;A PROCEDURE TO CREATE A BYTE ARRAY FOR MAKING A COMPOSITE MOVIE OF
+;MAGNETOCONVECTION OR DYNAMOS USING RESULTS FROM mm-files written by PROGRAM mag
+;
+;THIS PROCEDURE DISPLAYS COMPOSITE IMAGES OF Br ON THE OUTER SPHERICAL SURFACE
+;IN 3 MAP VIEWS: AITOFF, NORTH AND SOUTH POLAR, POSSIBLY ASSUMING
+;SOME LONGITUDINAL SYMMETRY. THE "NORTH" GEOMAGNETIC (DIPOLE) POLE IS ALSO PLOTTED
+;
+;THIS PROCEDURE ALSO CREATES JPG FILES OF THE IMAGES, FOR WEB MOVIES.
+;
+;NOTE: The main byte array containing the images has four dimensions in this
+;version, FRAMES=bytarr(nx,ny,3,nframes). This is for reading images from Mac-type
+;true-color screen displays, with image dimensions (nx,ny,3). The IDL keyword TRUE=3
+;then applies to the TVRD() and WRITE_JPEG commands. This may have to be modified
+;depending on the type of display the image is created on.
+;Consult IDL manual for details.
+;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+;SET UP X-WINDOW FOR COLOR DISPLAY
+SET_PLOT,'X'
+;SET DEVICE FOR COLOR on Mac
+DEVICE, RETAIN=2, DECOMPOSED=0
+
+;DEFINE THE VALUE OF PI:
+ PI= 3.1415926
+;ACCESS RGB-VALUES OF COLOR TABLE
+ COMMON COLORS,R_ORIG,G_ORIG,B_ORIG,R_CURR,G_CURR,B_CURR
+
+;THIS VERSION READS TIME AS A STRING VARIABLE
+TIME=' '
+
+;READS IN DATA FILENAME
+ FILENAME=''
+ GIFNAME=''
+ GIFFULLNAME=''
+ TSTEP=''
+ PRINT,'ENTER mm-type MOVIE DATA FILENAME:'
+ READ,FILENAME
+
+;OPEN DATA FILE AND READ PRELIMINARIES
+ OPENR, 1, FILENAME ;OPEN FILE FOR READ
+ RUNID=' '
+ READF,1,RUNID ;IDENTIFYING HEADER: RUNID
+ READF,1,TSCALE,VSCALE ;TIME SCALE, VELOCITY SCALE
+ READF,1,NR,NI,NPFULL,NSYM,NFRAME ;RECORD NUMBERS & FRAMES
+ PRINT,RUNID
+ PRINT,'NR,NI,NPFULL,NSYM,NFRAME= ',NR,NI,NPFULL,NSYM,NFRAME
+
+;DEFAULTS
+ NFSTOP=NFRAME-1
+ IINTER=0
+;INTERPOLATION AND FRAMES
+ PRINT,'Stop at frame NFSTOP ?'
+ READ,NFSTOP
+ PRINT,'Interpolate frames? No=0 Yes=1'
+ READ,IINTER
+;NUMBER OF FRAMES
+ NFR=(NFSTOP+2)*(IINTER+1)
+
+;NOTE: NSYM=N-FOLD SYMMETRY IN PHI:
+ NP=NPFULL/NSYM
+ NPFULL=NPFULL+1 ;WRAP-AROUND POINT ADDED
+
+;CHOOSE A VARIABLE
+ CHOICE=1 ;outer boundary Br
+;CHOOSE A DISPLAY
+ PRINT,'DISPLAY: Br at CMB Aitoff (top); N and S orthographic (bot)'
+ PRINT,'ENTER LONGITUDE (-180 to 180):'
+ READ, LONSHIFT & INCL=90
+ PRINT, 'WANT CONTINENTS? (yes=1):'
+ READ,ICONT
+ PRINT,'Select Polarity (Normal=1; Reversed=-1):'
+ READ,POLAR
+ IF POLAR NE -1 THEN POLAR=1
+
+;GIF-FILE OPTION
+ PRINT,'CREATE JPEG-FILES FOR WEB DISPLAY ? NO=0 YES=1'
+ READ,IGIF
+ IF IGIF EQ 1 THEN BEGIN
+ PRINT,'ENTER JPG-FILENAME (WITHOUT SUFFIXES):'
+ READ,GIFNAME
+ ENDIF
+
+;DECLARE COORDINATE ARRAYS
+ THET = FLTARR(NI)
+ DTHET = FLTARR(NPFULL,NI)
+ PHI = FLTARR(NPFULL)
+ LON = FLTARR(NPFULL)
+ LAT = FLTARR(NI)
+ COSPH2D = FLTARR(NPFULL,NI)
+ SINPH2D = FLTARR(NPFULL,NI)
+ COSTH2D = FLTARR(NPFULL,NI)
+ SINTH2D = FLTARR(NPFULL,NI)
+ DAREA = FLTARR(NPFULL,NI)
+
+;DEFINE DATA ARRAY SIZES
+ V1 = FLTARR(NP)
+ B1 = FLTARR(NP)
+ B0 = FLTARR(NP)
+ V = FLTARR(NPFULL,NI)
+ VI = FLTARR(NPFULL,NI)
+
+;DEFINE WINDOW SIZE
+ WXSIZE=300 & WYSIZE=400
+ PRINT,'SPECIFY (X,Y) WINDOW SIZE (eg, 300 400):'
+ READ,WXSIZE,WYSIZE
+
+;DIMENSION THE BYTE ARRAY FOR TRUE=3
+ FRAMES=BYTARR(WXSIZE,WYSIZE,3,NFR)
+
+;CREATE A WINDOW
+ WINDOW,2,XSIZE=WXSIZE,YSIZE=WYSIZE
+ WSET,2
+
+;SPECIFY LOCATIONS FOR COMBO PLOTS
+ XDIM=18.0 & YDIM=24.0
+ XC1=[1.5/XDIM,12./YDIM,17./XDIM,20./YDIM]
+ XC2=[1.5/XDIM,1.5/YDIM,17./XDIM,9.5/YDIM]
+ XY3=[ 0.5/XDIM, 2./YDIM, 8.5/XDIM,10./YDIM]
+ XT3=4.5/XDIM & YT3=10.2/YDIM
+ XY4=[ 9.5/XDIM, 2./YDIM,17.5/XDIM,10./YDIM]
+ XT4=13.5/XDIM & YT4=10.2/YDIM
+
+;READ IN COLATITUDES
+ READF,1,THET
+
+;create array PHI of longitudes of data points (in radians)
+ PI2NP = 2*PI/(NPFULL-1)
+ PHI=PI2NP*(FINDGEN(NPFULL))
+;converts PHI to array LON from -180 to 180 degrees
+ LON=180*PHI/PI - 180.0
+;converts THET into array LAT from 90 to -90 degrees
+ LAT=90-180*THET/PI
+;create 2D d(area) and trig arrays on unit sphere
+ FOR IT=0,NI-1 DO BEGIN
+ COSPH2D(*,IT)=COS(PHI(*))
+ SINPH2D(*,IT)=SIN(PHI(*))
+ ENDFOR
+ NF=NPFULL-1
+ FOR IP=0,NF DO BEGIN
+ COSTH2D(IP,*)=COS(THET)
+ SINTH2D(IP,*)=SIN(THET)
+ ENDFOR
+ FOR IP=0,NF DO BEGIN
+ DTHET(IP,*)=(SHIFT(THET,-1)-SHIFT(THET,1))/2
+ ENDFOR
+ DTHET(*,0)=DTHET(*,1) & DTHET(*,NI-1)=DTHET(*,NI-2)
+ DAREA=SINTH2D*DTHET*PI2NP
+;end of darea and trig array creation
+
+;SET UP CONTOUR LEVELS AND COLOR SCHEME
+
+;NUMBER OF CONTOURS
+ NLV=17
+;SCALE FACTOR
+ VSCALE=1.0
+;COLORS
+ CTABLE=39 ; 39: rainbow, 0: black/white
+ LOADCT,CTABLE ;LOADS COLOR TABLE NUMBER
+; SELF-DEFINED COLOR TABLE
+; BLACK COLOR AT START OF SPECTRUM
+ R_CURR(000)= 00 & G_CURR(000)= 00 & B_CURR(000)= 00
+; WHITE COLOR AT END OF SPECTRUM
+; R_CURR(255)=255 & G_CURR(255)=255 & B_CURR(255)=255
+; DARK GREY
+ R_CURR( 01)= 90 & G_CURR( 01)= 90 & B_CURR( 01)= 90
+; LIGHT GREY
+ R_CURR( 02)=190 & G_CURR( 02)=190 & B_CURR( 02)=190
+
+ R_CURR( 60)= 35 & G_CURR( 70)= 0 & B_CURR( 70)=230
+ R_CURR( 70)= 20 & G_CURR( 70)= 0 & B_CURR( 70)=245
+ R_CURR( 80)= 0 & G_CURR( 80)= 60 & B_CURR( 80)=255
+ R_CURR( 90)= 0 & G_CURR( 90)=110 & B_CURR( 90)=255
+ R_CURR(100)= 0 & G_CURR(100)=150 & B_CURR(100)=255
+ R_CURR(110)= 20 & G_CURR(110)=200 & B_CURR(110)=255
+ R_CURR(120)= 70 & G_CURR(120)=240 & B_CURR(120)=255
+ R_CURR(130)=150 & G_CURR(130)=255 & B_CURR(130)=240
+ R_CURR(140)=225 & G_CURR(140)=255 & B_CURR(140)=225
+
+ R_CURR(150)=255 & G_CURR(150)=255 & B_CURR(150)=185
+ R_CURR(160)=255 & G_CURR(160)=255 & B_CURR(160)= 75
+ R_CURR(170)=255 & G_CURR(170)=225 & B_CURR(170)= 0
+ R_CURR(180)=255 & G_CURR(180)=175 & B_CURR(180)= 0
+ R_CURR(190)=255 & G_CURR(190)=125 & B_CURR(190)= 0
+ R_CURR(192)=255 & G_CURR(192)= 75 & B_CURR(192)= 0
+ R_CURR(194)=255 & G_CURR(194)= 0 & B_CURR(194)= 0
+ R_CURR(196)=230 & G_CURR(196)= 0 & B_CURR(196)= 35
+
+ TVLCT,R_CURR,G_CURR,B_CURR
+
+ VCOLORS=[70,80,90,100,110,120,130,140,150,160,170,180,190,192,194,196,$
+ 196,70]
+
+; VCOLORS=[35,45,55,65,75,85,95,115,156,166,172,180,187,194,200,206,211]
+
+ !P.THICK=2 ;DRAWS THICK LINES
+ !P.NOERASE=1 ;PREVENT SPURIOUS ERASURE
+ !P.MULTI=[0,0,2] ;TWO PLOTS PER PAGE
+ !P.POSITION=[.025,.025,.975,.975] ;ENSURE SQUARE FIGURE
+ ICTC=0 ;NO TANGENT CYLINDER -- DEFAULT
+
+;************************************************************************
+;READ IN DATA ARRAYS FROM FORTRAN DATAFILE INTO IDL ARRAYS.
+;**************************************************
+;**************************************************
+;BEGIN TO MAKE MOVIE
+;**************************************************
+
+;MAIN FRAME LOOP
+ IFR=-1
+ IGIFFR=0
+FOR IFRAME=0,NFSTOP-1 DO BEGIN
+ IFR=IFR+1
+ READF,1,IDAT,TIME
+ PRINT,'FRAME= ',IDAT-1,' TIME= ',TIME
+
+;BUILD 2D ARRAYS
+
+;READ FIRST IMAGE
+ FOR IT=0,NI-1 DO BEGIN
+ READF,1,B0
+ FOR NS=0,NSYM-1 DO BEGIN
+ J0=NS*NP
+ J1=J0+NP-1
+ V(J0:J1,IT)=B0
+ ENDFOR
+ V(NPFULL-1,IT)=B0(0)
+ ENDFOR ;END OF FIRST IMAGE READ
+
+;READ SECOND IMAGE
+ FOR IT=0,NI-1 DO BEGIN
+ READF,1,B1
+ ENDFOR ;END OF SECOND IMAGE READ
+
+;READ THIRD IMAGE
+ FOR IT=0,NI-1 DO BEGIN
+ READF,1,B1
+ ENDFOR ;END OF THIRD IMAGE READ
+
+;AT FIRST FRAME, SET CONTOURING LEVELS
+LABEL0: IF IFRAME EQ 0 THEN BEGIN
+ VMAX=MAX(ABS(V)) / VSCALE
+ VSTEP=2*VMAX/(NLV-1)
+ VV = FINDGEN(NLV)*VSTEP-VMAX
+ VV(0) = VV(0) - 15*VSTEP ;DROP FIRST CONTOUR
+ ENDIF
+
+;begin geomagnetic pole calculation
+;spherical integration to get dipole moment components
+ MZ=TOTAL(V*COSTH2D*DAREA)-TOTAL(V(NF,*)*COSTH2D(NF,*)*DAREA(NF,*))
+ MX=TOTAL(V*SINTH2D*COSPH2D*DAREA) $
+ - TOTAL(V(NF,*)*SINTH2D(NF,*)*COSPH2D(NF,*)*DAREA(NF,*))
+ MY=TOTAL(V*COSTH2D*SINPH2D*DAREA) $
+ - TOTAL(V(NF,*)*COSTH2D(NF,*)*SINPH2D(NF,*)*DAREA(NF,*))
+ ME=SQRT(MX^2 + MY^2)
+;convert to geomagnetic components
+ MZ=-MZ & MX=-MX & MY=-MY
+;compute geomagnetic pole coordinates plat,plon
+ TILT=ATAN(ME,MZ) & PLAT=90-(180*TILT/PI)
+ TPHI=ATAN(MX,MY) & PLON=180*TPHI/PI
+;end of geomagnetic pole calculation
+
+
+;PLOT IMAGES
+
+;INTERPOLATED IMAGE
+IF IINTER EQ 1 AND IFRAME GT 0 THEN BEGIN
+ERASE
+;PLOT TOP IMAGE
+!P.POSITION=XC1
+TSTEP=STRTRIM(IFRAME,1)
+MAINID='CMB RADIAL FIELD AT t= '+ TSTEP
+MAP_SET,0,LONSHIFT,0,/AITOFF,TITLE=MAINID,/NOBORDER,/NOERASE
+CONTOUR,POLAR*V,LON,LAT,LEVELS=VV,/OVERPLOT,C_COLORS=VCOLORS,/CELL_FILL
+ IF (ICONT EQ 1) THEN MAP_CONTINENTS, MLINETHICK=1,COLOR=0
+ IF (ICONT NE 1) THEN MAP_GRID,COLOR=0,GLINESTYLE=1,GLINETHICK=1
+; PLOT TANGENT CYLINDER INTERSECTION
+ IF ICTC EQ 1 THEN BEGIN
+ TCLAT=69.6
+ PLOTS,2*FINDGEN(180),REPLICATE(TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ PLOTS,2*FINDGEN(180),REPLICATE(-TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ ENDIF
+; PLOT GEOMAGNETIC POLE POSITION
+ PLOTS,PLON,PLAT,PSYM=1,SYMSIZE=3,COLOR=250
+;PLOT LOWER IMAGES
+!P.POSITION=XY3
+MAP_SET,INCL,LONSHIFT,0,/ORTHOGRAPHIC,TITLE='NORTH',/NOERASE,/NOBORDER
+CONTOUR,POLAR*V,LON,LAT,LEVELS=VV,/OVERPLOT,C_COLORS=VCOLORS,/CELL_FILL
+ IF (ICONT EQ 1) THEN MAP_CONTINENTS, MLINETHICK=1,COLOR=0
+ IF (ICONT NE 1) THEN MAP_GRID,COLOR=0,GLINESTYLE=1,GLINETHICK=1
+; PLOT TANGENT CYLINDER INTERSECTION
+ IF ICTC EQ 1 THEN BEGIN
+ TCLAT=69.6
+ PLOTS,2*FINDGEN(180),REPLICATE(TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ PLOTS,2*FINDGEN(180),REPLICATE(-TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ ENDIF
+; PLOT GEOMAGNETIC POLE POSITION
+ PLOTS,PLON,PLAT,PSYM=1,SYMSIZE=3,COLOR=250
+!P.POSITION=XY4
+MAP_SET,-1*INCL,LONSHIFT,0,/ORTHOGRAPHIC,TITLE='SOUTH',/NOERASE,/NOBORDER
+CONTOUR,POLAR*V,LON,LAT,LEVELS=VV,/OVERPLOT,C_COLORS=VCOLORS,/CELL_FILL
+ IF (ICONT EQ 1) THEN MAP_CONTINENTS, MLINETHICK=1,COLOR=0
+ IF (ICONT NE 1) THEN MAP_GRID,COLOR=0,GLINESTYLE=1,GLINETHICK=1
+; PLOT TANGENT CYLINDER INTERSECTION
+ IF ICTC EQ 1 THEN BEGIN
+ TCLAT=69.6
+ PLOTS,2*FINDGEN(180),REPLICATE(TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ PLOTS,2*FINDGEN(180),REPLICATE(-TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ ENDIF
+; PLOT GEOMAGNETIC POLE POSITION
+ PLOTS,PLON,PLAT,PSYM=1,SYMSIZE=3,COLOR=250
+
+;CREATE FRAME
+ FRAMES(*,*,*,IFR)=TVRD(TRUE=3)
+ IFR=IFR+1
+
+;INCREMENT JPG FILE COUNTER
+ IF IGIF EQ 1 THEN IGIFFR=IGIFFR+1
+;CREATE JPG FILE
+ IF IGIFFR GT 0 THEN BEGIN
+ IF IGIFFR LT 10 THEN $
+ GIFFULLNAME=GIFNAME+'00'+STRTRIM(FIX(IGIFFR),2)+'.JPG' $
+ ELSE IF IGIFFR LT 100 THEN $
+ GIFFULLNAME=GIFNAME+'0'+STRTRIM(FIX(IGIFFR),2)+'.JPG'$
+ ELSE GIFFULLNAME=GIFNAME+STRTRIM(FIX(IGIFFR),2)+'.JPG'
+ WRITE_JPEG,GIFFULLNAME,TVRD(TRUE=3),TRUE=3
+ IGIFFR=IGIFFR+1
+ ENDIF
+ ENDIF
+;ENDIF ;END INTERPOLATED IMAGE
+
+
+;NEW IMAGE
+ERASE
+;PLOT TOP IMAGE
+!P.POSITION=XC1
+TSTEP=STRTRIM(IFRAME,1)
+MAINID='CMB RADIAL FIELD AT t= '+ TSTEP
+MAP_SET,0,LONSHIFT,0,/AITOFF,TITLE=MAINID,/NOBORDER,/NOERASE
+CONTOUR,POLAR*V,LON,LAT,LEVELS=VV,/OVERPLOT,C_COLORS=VCOLORS,/CELL_FILL
+ IF (ICONT EQ 1) THEN MAP_CONTINENTS, MLINETHICK=1,COLOR=0
+ IF (ICONT NE 1) THEN MAP_GRID,COLOR=0,GLINESTYLE=1,GLINETHICK=1
+; PLOT TANGENT CYLINDER INTERSECTION
+ IF ICTC EQ 1 THEN BEGIN
+ TCLAT=69.6
+ PLOTS,2*FINDGEN(180),REPLICATE(TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ PLOTS,2*FINDGEN(180),REPLICATE(-TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ ENDIF
+; PLOT GEOMAGNETIC POLE POSITION
+ PLOTS,PLON,PLAT,PSYM=1,SYMSIZE=3,COLOR=250
+;PLOT LOWER IMAGES
+!P.POSITION=XY3
+MAP_SET,INCL,LONSHIFT,0,/ORTHOGRAPHIC,TITLE='NORTH',/NOERAS,/NOBORDER
+CONTOUR,POLAR*V,LON,LAT,LEVELS=VV,/OVERPLOT,C_COLORS=VCOLORS,/CELL_FILL
+ IF (ICONT EQ 1) THEN MAP_CONTINENTS, MLINETHICK=1,COLOR=0
+ IF (ICONT NE 1) THEN MAP_GRID,COLOR=0,GLINESTYLE=1,GLINETHICK=1
+; PLOT TANGENT CYLINDER INTERSECTION
+ IF ICTC EQ 1 THEN BEGIN
+ TCLAT=69.6
+ PLOTS,2*FINDGEN(180),REPLICATE(TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ PLOTS,2*FINDGEN(180),REPLICATE(-TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ ENDIF
+; PLOT GEOMAGNETIC POLE POSITION
+ PLOTS,PLON,PLAT,PSYM=1,SYMSIZE=3,COLOR=250
+!P.POSITION=XY4
+MAP_SET,-1*INCL,LONSHIFT,0,/ORTHOGRAPHIC,TITLE='SOUTH',/NOERASE,/NOBORDER
+CONTOUR,POLAR*V,LON,LAT,LEVELS=VV,/OVERPLOT,C_COLORS=VCOLORS,/CELL_FILL
+ IF (ICONT EQ 1) THEN MAP_CONTINENTS, MLINETHICK=1,COLOR=0
+ IF (ICONT NE 1) THEN MAP_GRID,COLOR=0,GLINESTYLE=1,GLINETHICK=1
+; PLOT TANGENT CYLINDER INTERSECTION
+ IF ICTC EQ 1 THEN BEGIN
+ TCLAT=69.6
+ PLOTS,2*FINDGEN(180),REPLICATE(TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ PLOTS,2*FINDGEN(180),REPLICATE(-TCLAT,180),PSYM=3,SYMSIZE=.5,COLOR=0
+ ENDIF
+; PLOT GEOMAGNETIC POLE POSITION
+ PLOTS,PLON,PLAT,PSYM=1,SYMSIZE=3,COLOR=250
+;END NEW IMAGE
+
+;AT FIRST FRAME OPTION TO CHANGE SCALING OR PROJECTIONS
+ IF IFRAME EQ 0 THEN BEGIN
+ PRINT,'CHANGE SCALING OR ORIENTATION ? YES=1'
+ READ,IOPT
+ IF IOPT EQ 1 THEN BEGIN
+ PRINT,'SCALE FACTOR, LONGITUDE ?'
+ READ,VSCALE,LONSHIFT
+ PRINT, 'PLOT TANGENT CYLINDER INTERSECTION ? YES=1'
+ READ,ICTC
+ GOTO, LABEL0
+ END
+ END
+
+;LOAD IMAGE INTO BYTE ARRAY
+ FRAMES(*,*,*,IFR)=TVRD(TRUE=3)
+
+;INCREMENT JPG FILE COUNTER
+ IF IGIF EQ 1 THEN IGIFFR=IGIFFR+1
+;CREATE JPG FILE
+ IF IGIFFR GT 0 THEN BEGIN
+ IF IGIFFR LT 10 THEN $
+ GIFFULLNAME=GIFNAME+'00'+STRTRIM(FIX(IGIFFR),2)+'.JPG' $
+ ELSE IF IGIFFR LT 100 THEN $
+ GIFFULLNAME=GIFNAME+'0'+STRTRIM(FIX(IGIFFR),2)+'.JPG'$
+ ELSE GIFFULLNAME=GIFNAME+STRTRIM(FIX(IGIFFR),2)+'.JPG'
+ WRITE_JPEG,GIFFULLNAME,TVRD(TRUE=3),TRUE=3
+ ENDIF
+
+
+ENDFOR ;END OF MOVIE CREATION LOOP
+
+CLOSE,1 ;CLOSE DATA FILE
+
+
+
+;DISPLAY MOVIE
+
+LABEL20:
+
+PRINT,"MOVIE DISPLAY"
+
+REFRESH=0.1
+PRINT,'ENTER REFRESH TIME (nominal 0.1):'
+READ, REFRESH
+
+;RUN MOVIE
+FOR N=0,NFR-1 DO BEGIN
+ TV, FRAMES(*,*,*,N),TRUE=3
+ WAIT,REFRESH
+ENDFOR
+
+PRINT,'SHOW AGAIN? (YES=1):'
+READ,SHOWMORE
+IF (SHOWMORE EQ 1) THEN GOTO, LABEL20
+
+PRINT,"END MAGMOVCMB"
+!P.MULTI=0
+;END OF PROCEDURE magmovCMB
+
+END
Added: geodyn/3D/MAG/trunk/idl/magmovEQ3.pro
===================================================================
--- geodyn/3D/MAG/trunk/idl/magmovEQ3.pro 2007-04-24 19:40:51 UTC (rev 6655)
+++ geodyn/3D/MAG/trunk/idl/magmovEQ3.pro 2007-04-24 19:55:01 UTC (rev 6656)
@@ -0,0 +1,292 @@
+;A PROCEDURE TO CREATE A BYTE ARRAY FOR MAKING A MOVIE OF
+;MAGNETOCONVECTION OR DYNAMOS USING RESULTS FROM PROGRAM mag "me"
+;FILES, CONSISTING OF THREE SCALAR FIELDS IN THE EQUATORIAL PLANE
+;AS A FUNCTION OF TIME. THIS VERSION DISPLAYS ALL THREE IMAGES AT ONCE
+;
+;THIS PROCEDURE DISPLAYS RESULTS IN THE EQUATORIAL PLANE
+;POSSIBLY ASSUMING SOME LONGITUDINAL SYMMETRY.
+;
+;THIS PROCEDURE ALSO CREATES JPG FILES OF THE IMAGES, FOR WEB MOVIES.
+;
+;NOTE: The main byte array containing the images has four dimensions in this
+;version, FRAMES=bytarr(nx,ny,3,nframes). This is for reading images from Mac-type
+;true-color screen displays, with image dimensions (nx,ny,3). The IDL keyword TRUE=3
+;then applies to the TVRD() and WRITE_JPEG commands. This may have to be modified
+;depending on the type of display the image is created on.
+;Consult IDL manual for details.
+;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+;SET UP X-WINDOW FOR COLOR DISPLAY
+SET_PLOT,'X'
+;SET DEVICE FOR COLOR on Mac
+DEVICE, RETAIN=2, DECOMPOSED=0
+
+;DEFINE THE VALUE OF PI:
+ PI= 3.1415926
+;ACCESS RGB-VALUES OF COLOR TABLE
+ COMMON COLORS,R_ORIG,G_ORIG,B_ORIG,R_CURR,G_CURR,B_CURR
+
+;THIS VERSION READS TIME AS A STRING VARIABLE
+TIME=' '
+
+;READS IN DATA FILENAME
+ FILENAME=''
+ GIFNAME=''
+ GIFFULLNAME=''
+ TSTEP=''
+ PRINT,'ENTER MOVIE DATA FILENAME:'
+ READ,FILENAME
+
+;OPEN DATA FILE AND READ PRELIMINARIES
+ OPENR, 1, FILENAME ;OPEN FILE FOR READ
+ RUNID=' '
+ READF,1,RUNID ;IDENTIFYING HEADER: RUNID
+ READF,1,TSCALE,VSCALE ;TIME SCALE, VELOCITY SCALE
+ READF,1,NR,NI,NPFULL,NSYM,NFRAME ;RECORD NUMBERS & FRAMES
+ PRINT,RUNID
+ PRINT,'NR,NI,NPFULL,NSYM,NFRAME= ',NR,NI,NPFULL,NSYM,NFRAME
+
+;DEFAULTS
+ NFSTOP=NFRAME-1
+ IINTER=0
+;NUMBER OF FRAMES TO READ
+ PRINT,'Stop at frame NFSTOP ?'
+ READ,NFSTOP
+;NUMBER OF FRAMES
+ NFR=(NFSTOP+2)
+
+;NOTE: NSYM=N-FOLD SYMMETRY IN PHI:
+ NP=NPFULL/NSYM
+ NPFULL=NPFULL+1 ;WRAP-AROUND POINT ADDED
+
+;CHOOSE A DISPLAY
+ CHOICED=1 & LONSHIFT=0 ;equatorial plane
+ PRINT, "EQUATORIAL PLANE DISPLAYS"
+ POLAR=1
+
+
+;GIF-FILE OPTION
+ PRINT,'CREATE JPEG-FILES FOR WEB DISPLAY ? NO=0 YES=1'
+ READ,IGIF
+ IF IGIF EQ 1 THEN BEGIN
+ PRINT,'ENTER JPG-FILENAME (WITHOUT SUFFIXES):'
+ READ,GIFNAME
+ ENDIF
+
+;DECLARE COORDINATE ARRAYS
+ RAD = FLTARR(NR)
+ PHI = FLTARR(NPFULL)
+ LON = FLTARR(NPFULL)
+;DEFINE DATA ARRAY SIZES
+ A1 = FLTARR(NP)
+ A2 = FLTARR(NP)
+ A3 = FLTARR(NP)
+ T = FLTARR(NR,NPFULL)
+ B = FLTARR(NR,NPFULL)
+ W = FLTARR(NR,NPFULL)
+
+;DEFINE WINDOW SIZE
+ WXSIZE=600 & WYSIZE=200
+ PRINT,'SPECIFY WINDOW SIZE (eg, 600 200):'
+ READ,WXSIZE, WYSIZE
+
+;DIMENSION THE BYTE ARRAY FOR TRUE=3
+ FRAMES=BYTARR(WXSIZE,WYSIZE,3,NFR)
+
+;CREATE A WINDOW
+ WINDOW,2,XSIZE=WXSIZE,YSIZE=WYSIZE
+ WSET,2
+
+;READ IN RADII
+ READF,1,RAD
+
+;create array PHI of longitudes of data points (in radians)
+ PI2NP = 2*PI/(NPFULL-1)
+ PHI=PI2NP*(FINDGEN(NPFULL))
+;converts PHI to array LON from -180 to 180 degrees
+ LON=180*PHI/PI - 180.0
+
+;SET UP CONTOUR LEVELS AND COLOR SCHEME
+
+;NUMBER OF CONTOURS
+ NLV=17
+;SCALE FACTORS
+ TSCALE=1.0
+ BSCALE=1.0
+ WSCALE=1.0
+;COLORS
+ CTABLE=39 ; 39: rainbow, 0: black/white
+ LOADCT,CTABLE ;LOADS COLOR TABLE NUMBER
+; SELF-DEFINED COLOR TABLE
+; BLACK COLOR AT START OF SPECTRUM
+ R_CURR(000)= 00 & G_CURR(000)= 00 & B_CURR(000)= 00
+; WHITE COLOR AT END OF SPECTRUM
+; R_CURR(255)=255 & G_CURR(255)=255 & B_CURR(255)=255
+; DARK GREY
+ R_CURR( 01)= 90 & G_CURR( 01)= 90 & B_CURR( 01)= 90
+; LIGHT GREY
+ R_CURR( 02)=190 & G_CURR( 02)=190 & B_CURR( 02)=190
+
+ R_CURR( 60)= 35 & G_CURR( 60)= 0 & B_CURR( 60)=230
+ R_CURR( 70)= 20 & G_CURR( 70)= 0 & B_CURR( 70)=245
+ R_CURR( 80)= 0 & G_CURR( 80)= 60 & B_CURR( 80)=255
+ R_CURR( 90)= 0 & G_CURR( 90)=110 & B_CURR( 90)=255
+ R_CURR(100)= 0 & G_CURR(100)=150 & B_CURR(100)=255
+ R_CURR(110)= 20 & G_CURR(110)=200 & B_CURR(110)=255
+ R_CURR(120)= 70 & G_CURR(120)=240 & B_CURR(120)=255
+ R_CURR(130)=150 & G_CURR(130)=255 & B_CURR(130)=240
+ R_CURR(140)=225 & G_CURR(140)=255 & B_CURR(140)=225
+
+ R_CURR(150)=255 & G_CURR(150)=255 & B_CURR(150)=185
+ R_CURR(160)=255 & G_CURR(160)=255 & B_CURR(160)= 75
+ R_CURR(170)=255 & G_CURR(170)=225 & B_CURR(170)= 0
+ R_CURR(180)=255 & G_CURR(180)=175 & B_CURR(180)= 0
+ R_CURR(190)=255 & G_CURR(190)=125 & B_CURR(190)= 0
+ R_CURR(192)=255 & G_CURR(192)= 75 & B_CURR(192)= 0
+ R_CURR(194)=255 & G_CURR(194)= 0 & B_CURR(194)= 0
+ R_CURR(196)=230 & G_CURR(196)= 0 & B_CURR(196)= 35
+
+ TVLCT,R_CURR,G_CURR,B_CURR
+
+ CCOLORS=[70,80,90,100,110,120,130,140,150,160,170,180,190,192,194,196,$
+ 196,70]
+
+ !P.THICK=2 ;DRAWS THICK LINES
+ !P.MULTI=[0,3,0] ;THREE PLOTS IN A ROW
+; !P.POSITION=[.025,.025,.975,.975] ;ENSURE SQUARE FIGURE
+; XDIM=18.0 & YDIM=6.0
+; XC1=[1.5/XDIM,1.5/YDIM,16./XDIM,16./YDIM]
+
+;************************************************************************
+;READ IN DATA ARRAYS FROM MAG DATAFILE INTO IDL ARRAYS.
+;**************************************************
+;**************************************************
+;BEGIN TO MAKE MOVIE
+;**************************************************
+
+;MAIN FRAME LOOP
+ IFR=-1
+ IGIFFR=0
+FOR IFRAME=0,NFSTOP-1 DO BEGIN
+ IFR=IFR+1
+ READF,1,IDAT,TIME
+ PRINT,'FRAME= ',IDAT-1,' TIME= ',TIME
+
+;READ AND MAKE AN IMAGE
+ FOR IR=0,NR-1 DO BEGIN
+ READF,1,A1,A2,A3
+
+;ADD TO 2D ARRAYS
+ FOR NS=0,NSYM-1 DO BEGIN
+ J0=NS*NP
+ J1=J0+NP-1
+ T(IR,J0:J1)=A1
+ T(IR,NPFULL-1)=A1(0)
+ B(IR,J0:J1)=A2
+ B(IR,NPFULL-1)=A2(0)
+ W(IR,J0:J1)=A3
+ W(IR,NPFULL-1)=A3(0)
+ ENDFOR
+ ENDFOR ;END OF IMAGE BUILD
+
+;AT FIRST FRAME, SET CONTOURING LEVELS
+LABEL0: IF IFRAME EQ 0 THEN BEGIN
+ TMAX=MAX(ABS(T)) / TSCALE
+ TSTEP=2*TMAX/(NLV-1)
+ TT = FINDGEN(NLV)*TSTEP-TMAX
+ TT(0) = TT(0) - 10*TSTEP ;DROP FIRST CONTOUR
+ BMAX=MAX(ABS(B)) / BSCALE
+ BSTEP=2*BMAX/(NLV-1)
+ BB = FINDGEN(NLV)*BSTEP-BMAX
+ BB(0) = BB(0) - 10*BSTEP ;DROP FIRST CONTOUR
+ WMAX=MAX(ABS(W)) / WSCALE
+ WSTEP=2*WMAX/(NLV-1)
+ WW = FINDGEN(NLV)*WSTEP-WMAX
+ WW(0) = WW(0) - 10*WSTEP ;DROP FIRST CONTOUR
+ ENDIF
+
+; CONTOUR PLOTS
+ ERASE
+ MAINID='EQUATORIAL PLANE: TEMPERATURE'
+ POLAR_CONTOUR,POLAR*TRANSPOSE(T),PHI,RAD,LEVELS=TT,C_COLORS=CCOLORS,$
+ /FILL,XSTYLE=4,YSTYLE=4,TITLE=MAINID,CHARSIZE=1.5
+ POLYFILL,rad(nr-1)*cos(phi),rad(nr-1)*sin(phi),COLOR=0
+ OPLOT,/POLAR,REPLICATE(RAD(0),NPFULL),PHI,COLOR=255 ;DRAW OUTER EQUATOR
+ OPLOT,/POLAR,REPLICATE(RAD(NR-1),NPFULL),PHI,COLOR=255 ;DRAW OUTER EQUATOR
+
+ MAINID='AXIAL VORTICITY'
+ POLAR_CONTOUR,POLAR*TRANSPOSE(W),PHI,RAD,LEVELS=WW,C_COLORS=CCOLORS,$
+ /FILL,XSTYLE=4,YSTYLE=4,TITLE=MAINID,CHARSIZE=1.5
+ POLYFILL,rad(nr-1)*cos(phi),rad(nr-1)*sin(phi),COLOR=0
+ OPLOT,/POLAR,REPLICATE(RAD(0),NPFULL),PHI,COLOR=255 ;DRAW OUTER EQUATOR
+ OPLOT,/POLAR,REPLICATE(RAD(NR-1),NPFULL),PHI,COLOR=255 ;DRAW OUTER EQUATOR
+
+ TMSTEP=STRTRIM(IFRAME,1)
+ MAINT='AXIAL FIELD AT t= '+ TMSTEP
+ POLAR_CONTOUR,POLAR*TRANSPOSE(B),PHI,RAD,LEVELS=BB,C_COLORS=CCOLORS,$
+ /FILL,XSTYLE=4,YSTYLE=4,TITLE=MAINT,CHARSIZE=1.5
+ POLYFILL,rad(nr-1)*cos(phi),rad(nr-1)*sin(phi),COLOR=0
+ OPLOT,/POLAR,REPLICATE(RAD(0),NPFULL),PHI,COLOR=255 ;DRAW OUTER EQUATOR
+ OPLOT,/POLAR,REPLICATE(RAD(NR-1),NPFULL),PHI,COLOR=255 ;DRAW OUTER EQUATOR
+
+;AT FIRST FRAME OPTION TO CHANGE SCALING OR PROJECTIONS
+ IF IFRAME EQ 0 THEN BEGIN
+ PRINT,'CHANGE SCALING OR ORIENTATION ? YES=1'
+ READ,IOPT
+ IF IOPT EQ 1 THEN BEGIN
+ PRINT,'ENTER (Temp, Vort, B)-SCALE FACTORS, LONGITUDE ?'
+ READ,TSCALE,WSCALE,BSCALE,LONSHIFT
+ PRINT,'Select Polarity (Reversed=-1; others unchanged):'
+ READ,POLAR1 & IF POLAR1 EQ -1 THEN POLAR=POLAR1
+ GOTO, LABEL0
+ END
+ END
+
+;LOAD IMAGE INTO BYTE ARRAY
+ FRAMES(*,*,*,IFR)=TVRD(TRUE=3)
+
+;INCREMENT JPG FILE COUNTER
+ IF IGIF EQ 1 THEN IGIFFR=IGIFFR+1
+;CREATE JPG FILE
+ IF IGIFFR GT 0 THEN BEGIN
+ IF IGIFFR LT 10 THEN $
+ GIFFULLNAME=GIFNAME+'00'+STRTRIM(FIX(IGIFFR),2)+'.JPG' $
+ ELSE IF IGIFFR LT 100 THEN $
+ GIFFULLNAME=GIFNAME+'0'+STRTRIM(FIX(IGIFFR),2)+'.JPG'$
+ ELSE GIFFULLNAME=GIFNAME+STRTRIM(FIX(IGIFFR),2)+'.JPG'
+ WRITE_JPEG,GIFFULLNAME,TVRD(TRUE=3),TRUE=3
+ ENDIF
+
+
+ENDFOR ;END OF MOVIE CREATION LOOP
+
+CLOSE,1 ;CLOSE DATA FILE
+
+
+
+;DISPLAY MOVIE
+
+LABEL20:
+
+PRINT,"MOVIE DISPLAY"
+
+REFRESH=0.1
+PRINT,'ENTER REFRESH TIME (nominal 0.1):'
+READ, REFRESH
+
+;RUN MOVIE
+FOR N=0,NFR-1 DO BEGIN
+ TV, FRAMES(*,*,*,N),TRUE=3
+ WAIT,REFRESH
+ENDFOR
+
+PRINT,'SHOW AGAIN? (YES=1):'
+READ,SHOWMORE
+IF (SHOWMORE EQ 1) THEN GOTO, LABEL20
+
+PRINT,"END MAGMOVIE"
+
+;END OF PROCEDURE magmvovEQ3
+
+END
More information about the cig-commits
mailing list