[cig-commits] r6122 - geodyn/3D/MAG/trunk/idl
wei at geodynamics.org
wei at geodynamics.org
Tue Feb 27 11:52:56 PST 2007
Author: wei
Date: 2007-02-27 11:52:56 -0800 (Tue, 27 Feb 2007)
New Revision: 6122
Added:
geodyn/3D/MAG/trunk/idl/magline.pro
Log:
New IDL program to display data g-file
Added: geodyn/3D/MAG/trunk/idl/magline.pro
===================================================================
--- geodyn/3D/MAG/trunk/idl/magline.pro 2007-02-27 04:48:24 UTC (rev 6121)
+++ geodyn/3D/MAG/trunk/idl/magline.pro 2007-02-27 19:52:56 UTC (rev 6122)
@@ -0,0 +1,449 @@
+;PROCEDURE MAGLINE.PRO
+
+;A PROCEDURE TO DISPLAY RESULTS FROM MAG CALCULATIONS
+;IN THE FORM OF MAGNETIC FIELD LINES. THIS ROUTINE READS MAG
+; g-FILES AND PLOTS VARIOUS FIELD LINE IMAGES
+
+;SET UP X-WINDOW FOR COLOR DISPLAY
+SET_PLOT,'X'
+;SET DEVICE FOR COLOR on Mac
+DEVICE, RETAIN=2, DECOMPOSED=0
+;SET COLOR TABLE
+LOADCT,39
+
+;DEFINE THE VALUE OF PI:
+ PI= 3.1415926
+
+;READS IN DATA FILENAME
+ INFILE=''
+ PRINT,'ENTER DATA FILENAME:'
+ READ,INFILE
+
+;OPEN DATA FILE AND READ PRELIMINARIES
+ OPENR, 1,INFILE ;OPEN FILE FOR READ
+ RUNID=' '
+ READF,1,RUNID ;IDENTIFYING HEADER: RUNID
+ READF,1,TIME,NN,NI,NJ,NGRAD,NGCOLAT,NGLON,NSYM ;GRID PARAMETERS
+ READF,1,RA,EK,PR,PM,RADRATIO ;CONTROL PARAMETERS
+ PRINT,RUNID
+ RUNID=STRTRIM(RUNID)
+
+;NOTE: NSYM=N-FOLD SYMMETRY IN PHI
+;DEFINES ARRAY SIZES
+ NR=(NN+NGRAD-1)/NGRAD ;# RADIAL GRID POINTS
+ NT=(NI+NGCOLAT-1)/NGCOLAT ;# THETA GRID POINTS
+ NPFULL=(NJ+NGLON-1)/NGLON ;# TOTAL PHI GRID POINTS
+ NP=NPFULL/NSYM ;# PHI-POINTS IN DATA FILE
+ DPHI = 2*PI/NPFULL
+
+;DECLARE COORDINATE ARRAYS
+ THETA = FLTARR(NT)
+ LAT = FLTARR(NT)
+ PHI = FLTARR(NPFULL+1)
+ LON = FLTARR(NPFULL+1)
+
+;read in array THETA of colatitudes of data points (radians):
+ READF,1,THETA
+
+
+;converts THETA to an array LAT of latitudes 90 to -90 degs North to South:
+ LAT=90-180*THETA/PI
+;create array PHI of longitudes of data points (in radians)
+ PHI=DPHI*FINDGEN(NPFULL+1)
+ PHIMIN=PHI(0) & PHIMAX=PHI(NPFULL)
+
+;converts PHI to array LON from -180 to 180 degrees
+ LON=180*PHI/PI - 180.0
+
+;DIMENSION STATEMENTS FOR 2D AND 3D ARRAYS
+ R=FLTARR(NR) ;array of radii to spherical 2D surfaces
+ T1 =FLTARR(NP,NT) ;temperature on spherical 2D surface
+ VR1=FLTARR(NP,NT) ;radial velocity on spherical 2D surface
+ VT1=FLTARR(NP,NT) ;theta velocity on spherical 2D surface
+ VP1=FLTARR(NP,NT) ;phi velocity on spherical 2D surface
+ BR1=FLTARR(NP,NT) ;radial field on spherical 2D surface
+ BT1=FLTARR(NP,NT) ;theta field on spherical 2D surface
+ BP1=FLTARR(NP,NT) ;phi field on spherical 2D surface
+
+
+ BR=FLTARR(NPFULL+1,NT,NR) ;radial field in 3D shell
+ BT=FLTARR(NPFULL+1,NT,NR) ;theta field in 3D shell
+ BP=FLTARR(NPFULL+1,NT,NR) ;phi field in 3D shell
+ BAMP=FLTARR(NPFULL+1,NT,NR) ;field intenstiy in 3D shell
+
+
+ NLN=1000 ;field line dimension
+ X=FLTARR(NLN) ;field line x-position
+ Y=FLTARR(NLN) ;field line y-position
+ Z=FLTARR(NLN) ;field line z-position
+
+
+ XTH=FLTARR(NPFULL+1) ;core boundary coordinate
+ YTH=FLTARR(NPFULL+1) ;core boundary coordinate
+ ZTH=FLTARR(NPFULL+1) ;core boundary coordinate
+
+ XP=FLTARR(NT) ;core boundary coordinate
+ YP=FLTARR(NT) ;core boundary coordinate
+ ZP=FLTARR(NT) ;core boundary coordinate
+
+
+;************************************************************************
+;LOOP READS IN DATA ARRAYS FROM FORTRAN DATAFILE INTO IDL ARRAYS.
+;************************************************************************
+ FOR IR=0,NR-1 DO BEGIN
+
+ READF,1, format='(i4,e13.5)', KC, R1
+ PRINT, KC, R1 ;just for figuring out loops, remove later.
+ R(IR)=R1
+ READF,1,T1
+ READF,1,VR1
+ READF,1,VT1
+ READF,1,VP1
+ READF,1,BR1
+ READF,1,BT1
+ READF,1,BP1
+
+
+; ;MAKES 3D SHELL ARRAYS FROM
+; ;2D SPHERICAL SURFACE ARRAYS:
+
+ FOR J=0,NP-1 DO BEGIN
+ FOR NS=1,NSYM DO BEGIN
+ JNS=(NS-1)*NP+J
+ BR(JNS,*,IR)= BR1(J,*)
+ BT(JNS,*,IR)= BT1(J,*)
+ BP(JNS,*,IR)= BP1(J,*)
+ ENDFOR
+ ENDFOR
+ ENDFOR
+
+;INCLUDE ARRAY VALUES AT PHI=2*PI
+ BR(NPFULL,*,*)=BR(0,*,*)
+ BT(NPFULL,*,*)=BT(0,*,*)
+ BP(NPFULL,*,*)=BP(0,*,*)
+
+;***********************************************************************
+;END OF LOOP
+;***********************************************************************
+
+CLOSE,1 ;CLOSE DATA FILE
+
+
+;NON-DIMENSIONALIZES RADII VALUES SUCH THAT R_CMB = 1.0
+ RADTOP=R(0)
+ R=R/RADTOP
+
+;DO SOME SMOOTHING
+ PRINT,'ENTER SMOOTHING FACTOR, ODD (NOMINAL=3):'
+ READ,SMFACT
+ IF SMFACT GT 2 THEN BEGIN
+ BR=SMOOTH(BR,SMFACT)
+ BT=SMOOTH(BT,SMFACT)
+ BP=SMOOTH(BP,SMFACT)
+ ENDIF
+
+;FORM FIELD INTENSITY
+ BAMP=SQRT(BR*BR + BT*BT + BP*BP)
+
+;COLORS
+ CTABLE=39 ; 39: rainbow, 0: black/white
+ LOADCT,CTABLE ;LOADS COLOR TABLE NUMBER
+ ICOLOR=30 ;INNER CORE
+ OCOLOR=60 ;OUTER CORE
+ NCOLOR=220 ;NORMAL POLARITY
+ RCOLOR=220 ;REVERSE POLARITY
+ BACKCOLOR=0 ;BACKGROUND COLOR
+
+
+;SPHERICAL GRID SPACING FACTOR
+ DENGRID=2
+;CHARACTER SIZE
+; !Z.CHARSIZE=2.5
+ !X.CHARSIZE=2.5
+;FIELD LINE THICKNESS
+ LNTHICK=1
+;BACKGROUND
+ !P.BACKGROUND=BACKCOLOR
+;DEFAULT TITLE
+ PTITLE='' & PTITLE=RUNID
+
+;*************************************************************************
+;************** HERE THE PLOTTING MENU STARTS **************************
+;*************************************************************************
+
+
+LABEL0: PRINT,"ENTER OPTION:"
+ PRINT,"STOP=0; DRAW FIELD LINES=1; ERASE=2; COLORS=3;"
+ PRINT,"PRINT JPG=4; CHANGE TITLE=5"
+ READ,IOPTION
+
+ CASE IOPTION OF
+ 0: GOTO, LABEL99
+ 1: GOTO, LABEL1
+ 2: GOTO, LABEL2
+ 3: GOTO, LABEL3
+ 4: GOTO, LABEL4
+ 5: GOTO, LABEL5
+ ELSE: GOTO, LABEL0
+ ENDCASE
+
+LABEL2: ERASE
+ GOTO, LABEL0
+
+LABEL1: PRINT, "ENTER NUMBER OF FIELD LINES IN THETA, (eg, 10):"
+ READ, NLINET
+; PRINT, "ENTER NUMBER OF FIELD LINES IN LONGITUDE:"
+; READ, NLINEP
+ PRINT, "ENTER STARTING RADIUS (eg, 0.97):"
+ READ, RSTART
+ PRINT, "ENTER FIELDLINE STEP SIZE(eg, 0.05):"
+ READ, DS
+ PRINT, "ENTER X-AXIS ROTATION(eg, 90 or 30; 0 is distorted):"
+ READ, XAXISROT
+
+
+;OPEN THE Z-BUFFER
+ SET_PLOT,'Z'
+ ERASE
+ DEVICE,SET_RESOLUTION=[800,800]
+;FILL BACKGROUND
+ POLYFILL,[0,1,1,0],[0,0,1,1],/NORMAL,COLOR=BACKCOLOR
+
+;DRAW INNER CORE
+;DRAW INNER CORE LATITUDE CIRCLES; LOOP OVER THETA
+ FOR IT=0,NT-1,DENGRID DO BEGIN
+ XTH=0.36*SIN(THETA(IT))*COS(PHI)
+ YTH=0.36*SIN(THETA(IT))*SIN(PHI)
+ ZTH=REPLICATE(0.36*COS(THETA(IT)),NPFULL+1)
+;PLOT LATITUDE CIRCLE
+ PLOT_3DBOX,XTH,YTH,ZTH,/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=ICOLOR,GRIDSTYLE=1
+ ENDFOR
+
+;DRAW INNER CORE LONGITUDE CIRCLES; LOOP OVER PHI
+ FOR IP=0,NPFULL,DENGRID DO BEGIN
+ XP=0.36*SIN(THETA)*COS(PHI(IP))
+ YP=0.36*SIN(THETA)*SIN(PHI(IP))
+ ZP=0.36*COS(THETA)
+;PLOT LONGITUDE CIRCLE
+ PLOT_3DBOX,XP,YP,ZP,/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=ICOLOR,XTITLE=PTITLE,GRIDSTYLE=1
+ ENDFOR
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;MAIN DO LOOPS -- LOOP OVER STEPS IN THETA AND PHI
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+;USE TRIANGLULAR DISTRIBUTION OF FIELD LINES
+
+;SPECIFY STEP SIZE IN THETA
+;old TSTEP=PI/((NLINET+1)*2)
+ TSTEP=PI/(NLINET+1) ;new
+
+;LOOP IN THETA
+ FOR NLT=1,NLINET DO BEGIN
+
+;SPECIFY STEP SIZE IN PHI -- TRIANGULAR
+ NLINEP=4*NLINET*SIN(2*NLT/NLINET) ;new
+ PHISTEP=2*PI/(NLINEP)
+
+;LOOP IN PHI
+ FOR NLP=1,NLINEP DO BEGIN
+
+; INITIALIZE THE POSITIONS
+ RLN=RSTART
+ THETALN=NLT*TSTEP
+ PHILN=(NLP-1)*PHISTEP
+
+;FIND STARTING GRID POINT INDICES AND FIELD VALUES
+ IP=ROUND(PHILN/DPHI)
+ ITSORT=SORT(ABS(THETALN-THETA)) & IT=ITSORT(0)
+ IRSORT=SORT(ABS(RLN-R)) & IR=IRSORT(0)
+ BPLN=BP(IP,IT,IR)/BAMP(IP,IT,IR)
+ BTLN=BT(IP,IT,IR)/BAMP(IP,IT,IR)
+ BRLN=BR(IP,IT,IR)/BAMP(IP,IT,IR)
+
+;PRINT, "STARTING GRID NEIGHBORS (IP,IT,IR)
+; PRINT, IP,IT,IR
+
+;DETERMINE BRSIGN, OF STARTING BR, TO DIRECT INTEGRATION
+ IF BRLN GE 0 THEN BRSIGN=-1
+ IF BRLN LT 0 THEN BRSIGN=1
+
+;FORM STARTING CARTESIAN POINTS
+ Z(0)=RLN*COS(THETALN)
+ X(0)=RLN*SIN(THETALN)*COS(PHILN)
+ Y(0)=RLN*SIN(THETALN)*SIN(PHILN)
+
+;SPECIFY FIELD LINE POLARITY (GEOMAGNETIC CONVENTION)
+ POLARITY=-Z(0)*BRSIGN & POLARITY=POLARITY/ABS(POLARITY)
+;SELECT COLOR
+ IF POLARITY GE 0 THEN FCOLOR=NCOLOR
+ IF POLARITY LT 0 THEN FCOLOR=RCOLOR
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;BEGIN DO LOOP OVER FIELD LINE POINTS
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+ FOR ILN=1,NLN-2 DO BEGIN
+
+;FIND NEAREST GRID POINT INDICES
+ IPSORT=SORT(ABS(PHI-PHILN)) & IP=IPSORT(0)
+ ITSORT=SORT(ABS(THETA-THETALN)) & IT=ITSORT(0)
+ IRSORT=SORT(ABS(R-RLN)) & IR=IRSORT(0)
+
+;FORM FIELD COMPONENTS
+ BPLN=BP(IP,IT,IR)/BAMP(IP,IT,IR)
+ BTLN=BT(IP,IT,IR)/BAMP(IP,IT,IR)
+ BRLN=BR(IP,IT,IR)/BAMP(IP,IT,IR)
+
+;ADVANCE FIELD LINE
+ PHILN=PHILN + BRSIGN*(BPLN*DS/(RLN*SIN(THETALN)))
+ THETALN=THETALN + BRSIGN*(BTLN*DS/(RLN))
+ RLN=RLN + BRSIGN*BRLN*DS
+
+;APPLY BOUND CHECKS
+ IF THETALN LT 0.01 THEN THETALN=0.01
+ IF THETALN GT PI-0.01 THEN THETALN=PI-0.01
+ IF PHILN GT PHIMAX THEN PHILN=PHILN-PHIMAX
+ IF PHILN LT 0 THEN PHILN=PHIMAX+PHILN
+
+;FORM CARTESIAN POINTS
+ Z(ILN)=RLN*COS(THETALN)
+ X(ILN)=RLN*SIN(THETALN)*COS(PHILN)
+ Y(ILN)=RLN*SIN(THETALN)*SIN(PHILN)
+
+;EXIT CHECKS
+ IF RLN GT 1 THEN GOTO, LABEL88
+ IF RLN LT 0.34 THEN GOTO, LABEL89
+
+;END OF FIELD LINE DO LOOP
+ ENDFOR
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; PLOT FIELD LINE
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+LABEL88:
+;PRINT, "FIELDLINE EXITS CORE AT ILN=",ILN
+; overplot field line
+ PLOT_3DBOX,X(0:ILN),Y(0:ILN),Z(0:ILN),/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=FCOLOR,THICK=LNTHICK,GRIDSTYLE=1
+
+GOTO, LABEL10
+
+LABEL89:
+;PRINT, "FIELDLINE ENTERS INNER CORE AT ILN=",ILN
+; overplot field line
+ PLOT_3DBOX,X(0:ILN),Y(0:ILN),Z(0:ILN),/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=FCOLOR,THICK=LNTHICK,GRIDSTYLE=1
+
+GOTO, LABEL10
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;END OF FIELD LINE INRECMENT LOOPS
+
+LABEL10:
+ ENDFOR
+ ENDFOR
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+
+;DRAW OUTER CORE
+
+;DRAW OUTER CORE LATITUDE CIRCLES; LOOP OVER THETA
+ FOR IT=0,4 DO BEGIN
+ THETAC=IT*PI/4
+ XTH=SIN(THETAC)*COS(PHI)
+ YTH=SIN(THETAC)*SIN(PHI)
+ ZTH=REPLICATE(COS(THETAC),NPFULL+1)
+;PLOT LATITUDE CIRCLE
+ PLOT_3DBOX,XTH,YTH,ZTH,/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=OCOLOR,GRIDSTYLE=1
+ ENDFOR
+
+;DRAW OUTER CORE LONGITUDE CIRCLES; LOOP OVER PHI
+ FOR IP=0,8 DO BEGIN
+ PHIC=IP*PI/4
+ XP=SIN(THETA)*COS(PHIC)
+ YP=SIN(THETA)*SIN(PHIC)
+ ZP=COS(THETA)
+;PLOT LONGITUDE CIRCLE
+ PLOT_3DBOX,XP,YP,ZP,/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=OCOLOR,GRIDSTYLE=1
+
+;OVERPLOT TO ELIMINATE GRID
+ PLOT_3DBOX,XP[0:1],YP[0:1],-ZP[0:1],/NOERASE,$
+ XRANGE=[-1,1],YRANGE=[-1,1],ZRANGE=[-1,1],$
+ AX=XAXISROT,COLOR=BACKCOLOR,GRIDSTYLE=1
+
+ ENDFOR
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+;READ IMAGE FROM Z-BUFFER
+ IMAGE=TVRD(TRUE=3)
+
+;SET DISPLAY TO SCREEN
+ SET_PLOT,'X'
+
+;SPECIFY WINDOW INFORMATION
+ XWIND=800
+ YWIND=800
+ WINDOW,0,XSIZE=XWIND,YSIZE=YWIND,TITLE='GRAPHIC WINDOW'
+
+;DISPLAY IMAGE
+ TVSCL,IMAGE
+
+;RETURN TO MAIN OPTIONS
+GOTO, LABEL0
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+
+LABEL3:
+PRINT, "ENTER COLOR TABLE (B&W=0; RAINBOW&WHITE=39):"
+READ,CTABLE & LOADCT,CTABLE
+PRINT, FORMAT='("CURRENT NORMAL,REVERSE,INNER CORE,OUTER CORE=",4I4)'$
+ ,NCOLOR,RCOLOR,ICOLOR,OCOLOR
+PRINT, "ENTER FOUR NEW COLORS:"
+READ, NCOLOR,RCOLOR,ICOLOR,OCOLOR
+PRINT,"ENTER BACKGROUND COLOR:"
+READ,BACKCOLOR & !P.BACKGROUND=BACKCOLOR
+PRINT, "CURRENT LINE THICKNESS=" & PRINT, LNTHICK
+PRINT, "ENTER NEW THICKNESS:"
+READ,LNTHICK
+PRINT, "CURRENT SPHERE GRID SPACING FACTOR=" & PRINT,DENGRID
+PRINT,"ENTER NEW FACTOR (INTEGER):"
+READ,DENGRID
+GOTO, LABEL0
+
+LABEL4:
+OUTFILE=''
+PRINT,"ENTER JPG FILENAME:"
+READ,OUTFILE
+;READS A 24-BIT SCREEN IMAGE, CONVERTS TO JPEG (IDL NO LONGER SUPPORTS GIF)
+AJPEG=TVRD(TRUE=3)
+WRITE_JPEG,OUTFILE,AJPEG,TRUE=3
+GOTO, LABEL0
+
+LABEL5:
+PRINT,"CURRENT TITLE:" & PRINT,PTITLE
+PRINT, "ENTER NEW TITLE:"
+READ,PTITLE
+GOTO, LABEL0
+
+
+LABEL99: PRINT, "EXIT FROM PROGRAM MAGLINE"
+
+END
+
More information about the cig-commits
mailing list