; FILE: butterfly.idl  [IDL source code]
;===========================================================
; LINEAR SPHERICAL MEAN-FIELD DYNAMO
;-----------------------------------------------------------
; APAS7500 Fall 1997  Problem IV.5.1 through 5
; see Section IV.5.2 of notes
;-----------------------------------------------------------
; This codes produces butterfly diagrams of B_r and
; B_phi at depths rbr and rbphi, as set by the user (see below).
; The solution is reconstructed from a solution
; produced by the fortran code sphdyn_eig.f
;-----------------------------------------------------------
; HARDCOPY: when producing a postscript hardcopy, before
; running the code (.run medplane.idl) make sure to set
;
; device,/portrait,yoffset=5.,xsize=18,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
 rbphi=0.7              ; depth of toroidal field cut
 rbr  =1.0              ; depth of radial field cut

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

 title0=string(' ',format='(a64)') & ctime=1.e0
 readu,1,title0 
 print,'Using Input file: ',solfile
 print,'Run title:        ',title0
 idum=lonarr(4)
 readu,1,idum
    numnp =idum(0)
    numel =idum(1)
    nelx  =idum(2)
    nely  =idum(3)
 readu,1,calpha,comega
 nt=50
 tout=2.
 nnodes=numnp/2

xx   =fltarr(nelx+1)
yy   =fltarr(nely+1)
b_r =fltarr(nt+1,nelx+1)
b_y =fltarr(nt+1,nelx+1)
aquad=fltarr(nelx+1)
daqdt=fltarr(nelx+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)
time=findgen(nt+1)/(nt)*tout

readu,1,xx
readu,1,yy
xxa=acos(xx)
;---------- find boundaries of layer and top
iblay=0
ibp=0
for k=1,nely-1 do if yy(k+1) gt rbphi  and yy(k) le rbphi  then ibp=k
for k=1,nely-1 do if yy(k+1) gt rbr    and yy(k) le rbr    then ibr=k
print, 'Bphi Butterfly diagram constructed at r/R=',yy(ibp)
print, 'Br   Butterfly diagram constructed at r/R=',yy(ibr)

;-----------------------------------------------------------------------
; 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]

;-----------------------------------------------------------------------
; plot grid

pos1 =[0.15,0.50,0.95,0.85]
pos2 =[0.15,0.10,0.95,0.45]

;-----------------------------------------------------------------------------
; prepare output file for movie

nnx=800
nny=800
;window,14,retain=2,xpos=50,ypos=200,xs=nnx,ys=nny

;-----------------------------------------------------------------------------
; loop over timesteps

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

;---------- Loop over phases

    for it=0,nt do begin

    phase=time(it)*2.*!PI
    sinphase=sin(phase)
    cosphase=cos(phase)

;---------- Reconstruct and extract toroidal field below interface
    for j=0,nelx do begin
       rsin=yy(ibp-1)*sqrt(1.-xx(j)^2)
       b_y(it,j)  =qr2(j,ibp-1)*cosphase-qi2(j,ibp-1)*sinphase
    endfor

;---------- Reconstruct vector potential
    for j=0,nelx do begin
       rsin=yy(ibr)*sqrt(1.-xx(j)^2)
       aquad(j)=rsin*(qr1(j,ibr)*cosphase-qi1(j,ibr)*sinphase)
    endfor

;---------- Compute vector potential angular derivative (to get B_r)
    for j=1,nelx-1 do begin
       sinthp=sqrt(1.-xx(j+1)^2)
       sinthm=sqrt(1.-xx(j-1)^2)
       delth=xxa(j+1)-xxa(j-1)
       daqdt(j)=(aquad(j+1)*sinthp-aquad(j-1)*sinthm)/delth
    endfor
    sinthp=sqrt(1.-xx(1)^2)
    sinthm=sqrt(1.-xx(0)^2)
    delth=xxa(1)-xxa(0)
    daqdt(0)=(aquad(1)*sinthp-aquad(0)*sinthm)/delth

;---------- Now finally calculate B_r at surface
    for j=1,nelx-1 do begin
       rsin=1.0*sqrt(1.-xx(j)^2)
       b_r(it,j)=daqdt(j)/rsin
    endfor
    thmid=0.5*(xxa(1)+xxa(0))
    rsinthmid=1.0*sin(thmid)
    b_r(it,0)=daqdt(0)/rsinthmid
    b_r(it,nelx)=-b_r(it,0)

endfor

;---------- Now normalize fields to [-1.,1.]
maxbr=max(b_r)
maxby=max(b_y)
b_r=b_r/maxbr
b_y=b_y/maxby

title1='B!9P!6(r/R='+string(yy(ibp),format='(f4.2)')+'), max(B!9P!6)='$
       +string(maxby,format='(e8.2)')
title2='Br(r/R='+string(yy(ibr),format='(f4.2)')+'), max(Br)='$
       +string(maxbr,format='(e8.2)')
; plot butterfly diagrams
;------------------------------------------------------------------------------
    !X.STYLE=1 & !X.THICK=2
    !Y.STYLE=1 & !Y.THICK=2
;---------- 1. plot color-filled contours for B-phi
    contour,pos=pos1,levels=lev,/noclip,$
       c_linestyle=clst,xcharsize=0.0000000001,charsize=1.25,$
       xrange=[0.0,tout],yrange=[-1.,1.],xticklen=0.04,$
       ytitle='!6(latitude)/(!7p!6/2)',$
       title='!6 ',$
       b_y,time,(2.*xxa/!PI-1.)
    xyouts,0.,1.5,title0
    xyouts,0.,1.1,title1
    relab='!6Re(!7k!n!6)='+string(reomega,format='(e9.2)')
    imlab='!6Im(!7k!n!6)='+string(imomega,format='(e9.2)')
    colab='!6C!7X!6='+string(comega,format='(e9.2)')
    calab='!6C!7a!6='+string(calpha,format='(e9.2)')
    xyouts,0.0,1.3,calab
    xyouts,0.5,1.3,colab
    xyouts,1.0,1.3,relab
    xyouts,1.5,1.3,imlab

    contour,/noerase,pos=pos2,levels=lev,/noclip,$
       c_linestyle=clst,charsize=1.25,$
       xrange=[0.0,tout],yrange=[-1.,1.],xticklen=0.04,$
       ytitle='!6(latitude)/(!7p!6/2)',$
       xtitle='!9P!6/(2!7p!6)',$
       b_r,time,(2.*xxa/!PI-1.)
    xyouts,0.,1.1,title2

end