; FILE: medplane.idl  [IDL source code]
;===========================================================
; LINEAR SPHERICAL DYNAMO
;-----------------------------------------------------------
; APAS7500 Fall 1997  Problem IV.5.1 through 5
; see Section 5.2 of notes
;-----------------------------------------------------------
; This codes reconstruct solution from the eigenvalue/vector
; produced by the fortran code sphdyn_eig.f
;
; Solution is plotted as B_phi isocontours and poloidal
; fiedlines in meridional quadrant, over half a dynamo
; cycle at phase intervals of pi/6
;-----------------------------------------------------------
; HARDCOPY: when producing a postscript hardcopy, before
; running the code (.run medplane.idl) make sure to set
;
; device,/landscape,xsize=24,ysize=18
;-----------------------------------------------------------
; INPUT: this code requires an input file, the name of which
; being set immediately below
;===========================================================
 solfile='solution.i3e' ; the name of your output file

;***** you shouldn't have to touch anything below this point *****
 !ORDER=1
 close,1 & openr,1,/f77_unformatted,solfile

 title0=string(' ',format='(a64)')

 readu,1,title0 & print,title0
 idum=lonarr(4)
 readu,1,idum
    numnp =idum(0) & numel =idum(1) & nelx  =idum(2) & nely  =idum(3)
 readu,1,calpha,comega

xx   =fltarr(nelx+1)
yy   =fltarr(nely+1)
b_y =fltarr(nelx+1,nely+1)
yc   =fltarr(nelx+1,nely+1)
xc   =fltarr(nelx+1,nely+1)
bquad=fltarr(nelx+1,nely+1)
qr1=fltarr(nelx+1,nely+1)
qi1=fltarr(nelx+1,nely+1)
qr2=fltarr(nelx+1,nely+1)
qi2=fltarr(nelx+1,nely+1)

readu,1,xx
readu,1,yy

for k=1,nely-1 do if yy(k+1) gt 1.0 and yy(k) le 1.0 then itop=k
;-----------------------------------------------------------------------
; select contour interval 
lev=2.*findgen(16)/15.-1.
clst=[1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0]

;--------- reconstruct grid in real [x,y] space
for i=0,nely do begin
   for j=0,nelx do begin
      sin1=sqrt(1.-xx(j)^2)
      xc(j,i)=yy(i)*sin1
      yc(j,i)=yy(i)*xx(j)
   endfor
endfor

pos1 =[0.1,0.65,0.225,0.90]
pos2 =[0.4,0.65,0.525,0.90]
pos3 =[0.7,0.65,0.825,0.90]
pos4 =[0.1,0.35,0.225,0.60]
pos5 =[0.4,0.35,0.525,0.60]
pos6 =[0.7,0.35,0.825,0.60]
pos7 =[0.1,0.05,0.225,0.30]

;---------- read results for eigenvector and eigenvalue
vdum=fltarr(3)
readu,1,vdum
   reomega=vdum(0) & imomega=vdum(1)
readu,1,qr1
readu,1,qi1
readu,1,qr2
readu,1,qi2

relab='!6Re(!7k!n!6)'
imlab='!6Im(!7k!n!6)'
colab='!6C!7X!6'
calab='!6C!7a!6'
relab2='!6= '+string(reomega,format='(f8.3)')
imlab2='!6= '+string(imomega,format='(f8.3)')
colab2='!6= '+string(comega,format='(f8.1)')
calab2='!6= '+string(calpha,format='(f8.1)')
;---------- Loop over phases

nt=7
time=1.*findgen(nt)/(nt-1)
for it=1,nt do begin

    if it eq 1 then begin
       posa=pos1
       phaselab='!6(A) !9P!6=0'
       erase
    endif
    if it eq 2 then begin
       posa=pos2
       phaselab='!6(B) !9P!6=!7p!6/6'
    endif
    if it eq 3 then begin
       posa=pos3
       phaselab='!6(C) !9P!6=!7p!6/3'
    endif
    if it eq 4 then begin
       posa=pos4
       phaselab='!6(D) !9P!6=!7p!6/2'
    endif
    if it eq 5 then begin
       posa=pos5
       phaselab='!6(E) !9P!6=2!7p!6/3'
    endif
    if it eq 6 then begin
       posa=pos6
       phaselab='!6(F) !9P!6=5!7p!6/6'
    endif
    if it eq 7 then begin
       posa=pos7
       phaselab='!6(G) !9P!6=!7p!6'
    endif

    phase=time(it-1)*!PI+!PI/6.
    sinphase=sin(phase)
    cosphase=cos(phase)

;---------- Reconstruct toroidal and poloidal fields
    for i=0,nely do begin
       for j=0,nelx do begin
          rsin=yy(i)*sqrt(1.-xx(j)^2)
          bquad(j,i)=rsin*(qr1(j,i)*cosphase-qi1(j,i)*sinphase)
          b_y(j,i)  =qr2(j,i)*cosphase-qi2(j,i)*sinphase
       endfor
    endfor

;---------- Now normalize fields to [-1.,1.]
    if it eq 1 then begin
       maxbq=max(bquad)
       maxby=max(b_y)
    endif
    bquad=bquad/maxbq
    b_y=b_y/maxby

;---------- now make the plots
    !X.STYLE=4 & !X.THICK=2 & !Y.STYLE=4 & !Y.THICK=2
;---------- 1. first isocontours of toroidal field strength
    contour,/noerase,/follow,pos=posa,levels=lev,/noclip,$
       c_labels=[0,0,0], c_linestyle=clst,$
       xrange=[-1.0,0.],yrange=[0.0,1.0],$
       b_y(nelx/2:nelx,0:itop),$
       -xc(nelx/2:nelx,0:itop),yc(nelx/2:nelx,0:itop)
   oplot,/noclip,-xc(nelx/2:nelx,itop:itop),yc(nelx/2:nelx,itop:itop)
   oplot,/noclip,[-1.,0.,0.],[0.,0.,1.]
;  if it eq 1 then xyouts,-1.,1.25,title0
   xyouts,-1.0,1.05,charsize=1.0,phaselab

;---------- 2. now poloidal fieldlines
   contour,/follow,/noclip,/noerase,pos=posa,levels=lev,$
      xrange=[-1.0,0.],yrange=[0.,1.0],$
       c_labels=[0,0,0],c_linestyle=clst,$ 
      bquad(nelx/2:nelx,0:itop+3),$
      xc(nelx/2:nelx,0:itop+3),yc(nelx/2:nelx,0:itop+3)
   oplot,/noclip,xc(nelx/2:nelx,itop:itop),yc(nelx/2:nelx,itop:itop),thick=2
   oplot,/noclip,[1.,0.],[0.,0.],thick=2

endfor
xyouts,charsize=1.25,1.5,0.8,colab
xyouts,charsize=1.25,1.5,0.6,calab
xyouts,charsize=1.25,1.5,0.4,relab
xyouts,charsize=1.25,1.5,0.2,imlab
xyouts,charsize=1.25,2.5,0.8,colab2
xyouts,charsize=1.25,2.5,0.6,calab2
xyouts,charsize=1.25,2.5,0.4,relab2
xyouts,charsize=1.25,2.5,0.2,imlab2
end