; FILE: helio14.idl  [IDL source code]
;===========================================================
; 2D AXISYMMETRIC MEAN-FIELD-LIKE DYNAMO
; 2014 HELIOPHYSICS SUMMER SCHOOL
;-----------------------------------------------------------
; These IDL procedures read in, analyses and displays various
; pre-calculated 2D axisymmetric kinematic mean-field dynamo
; solutions.
;-----------------------------------------------------------
; IDL PROCEDURES INCLUDED IN THIS FILE: 
; 
; mkcolort_pro.idl   : builds a custom false-color table (8-bit)
; colorscale_pro.idl : draws a color bar
; viewdynamo_pro.idl : the main analysis procedure
;-----------------------------------------------------------
; USAGE (see also section 1.6 in Lab Notes):
;
; first compile on IDL command line: 
;    IDL> .run helio14.idl
; then run on IDL command like:
;    IDL> .run viewdynamo
;
; This procedure REQUIRES a parameter file called param.txt, comprised
; of three text lines as follows:
;
;   solfile
;   rbphi
;   latbphi
;
;   where solfile is a character string giving the name of the dynamo
;   solution to be analyzed, rbphi is the radius at which time-latitude
;   diagram is to be constructed, and latbphi is the latitude at which
;   time-radius diagram is to be constructed
;===========================================================
PRO makecolortable
;--------------------------------------------
; construct convenient color table
;--------------------------------------------
loadct,0
rr=fltarr(32)
gg=fltarr(32)
bb=fltarr(32)
rr(20)=10 & gg(20)=10 & bb(20)=10
;---------- first made shade-o-yellow and shade-o-green scales
for k=1,6 do begin
   rr(20+k)=k*30 +40
   gg(20+k)=k*30 +40
   bb(20+k)=0 
   rr(20-k)=0
   gg(20-k)=k*30 +30
   bb(20-k)=0 
endfor
;---------- now add transition to red and cyan at upper ends
   rr(27)=240 & gg(27)=150 & bb(27)=0
   rr(28)=240 & gg(28)=100 & bb(28)=0
   rr(29)=240 & gg(29)=50  & bb(29)=0
   rr(30)=240 & gg(30)=0   & bb(30)=0
   rr(13)=0   & gg(13)=220 & bb(13)=120
   rr(12)=0   & gg(12)=170 & bb(12)=170
   rr(11)=0   & gg(11)=120 & bb(11)=220
   rr(10)=0   & gg(10)=0   & bb(10)=250
;---------- now make a few convenient additional colors
 rr(1)=240 & gg(1)=100 & bb(1)=0   ; orange
 rr(2)=150 & gg(2)=150 & bb(2)=150 ; gray
 rr(3)=150 & gg(3)=0   & bb(3)=250 ; purple
 rr(4)=0   & gg(4)=200 & bb(4)=200 ; cyan
 rr(5)=250 & gg(5)=250 & bb(5)=0   ; yellow
 rr(6)=0   & gg(6)=200 & bb(6)=0   ; green
 rr(7)=250 & gg(7)=0   & bb(7)=0   ; red
 rr(8)=0   & gg(8)=0   & bb(8)=250 ; blue
 rr(9)=255 & gg(9)=255 & bb(9)=255 ; white

tvlct,rr,gg,bb

return
end

;===========================================================
PRO colorscale,pos1,lev,cin,xi,xo,yi,yo,thk

;   prepare matrix for color scale
 eps=0.001
 cscal=intarr(4,24)
 xscal=fltarr(4,24)
 xscal(0:0,*)=xi
 xscal(1:1,*)=xi+eps
 xscal(2:2,*)=xo-eps
 xscal(3:3,*)=xo
 yscal=fltarr(4,24)
 for i=0,23 do yscal(*,i:i)=yi+(yo-yi)*i/23.
 yscal(*,0:0)=yi
 yscal(*,1:1)=yi+eps
 yscal(*,22:22)=yo-eps
 yscal(*,23:23)=yo
 levelc=cin
 cin2=intarr(21)
 cin2(0:20)=cin(1:21)-1
 cin2(0)=10 
 for j=1,22 do begin
    cscal(1:2,j)=j+9
    cscal(0:0,j)=-1
    cscal(3:3,j)=-1
 endfor
 cscal(0:3,0:0)=-1.e2
 cscal(0:3,23:23)=-1.e2

; flip colorscale
   !X.STYLE=5 & !Y.STYLE=5
;   3. color scale
   contour,/follow,/noclip,pos=pos1,levels=levelc,/noerase,$
      /fill,c_color=cin2,$
      xrange=[0.,1.],yrange=[0.,1.],$
      cscal,xscal,yscal
   ymid0=0.5*(yi+yo)
   eps=(yo-yi)/20.
;  lmin=string(alog10(abs(lev(0))),format='(i2)')
;  lmax=string(alog10(abs(lev(21))),format='(i2)')
;  lmid=string(alog10(lev(11)),format='(i2)')
;  xyouts,xo*1.01,yi+eps,charthick=2,charsize=1.,'!3-'+lmin
;  xyouts,xo*1.01,yo-eps,charthick=2,charsize=1.,'!3+'+lmax
;  xyouts,xo*1.01,ymid0,charthick=2,charsize=1.,'!9+!3'+lmid
   oplot,/noclip,thick=thk,[xi,xi,xo,xo,xi],[yi,yo,yo,yi,yi]
   !X.STYLE=1 & !Y.STYLE=1

;------------------------------------------------------------------------------
return
end

;===========================================================
PRO viewdynamo

; read in input parameters
 close,2 & openr,2,'param.txt'
 solfile=string(' ',format='(a30)') & ctime=1.e0
 readf,2,solfile
 readf,2,rbphi
 readf,2,latbphi
 close,2
 print,'Using Input file: ',solfile

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

;---------- read and print file header
 title0=string(' ',format='(a80)') & ctime=1.e0
 readu,1,title0 
;print,'Run title:        ',title0
 idum=lonarr(6)
 vdum=fltarr(8)
 readu,1,idum
    numnp =idum(0)
    numel =idum(1)
    nelx  =idum(2)
    nely  =idum(3)
    nt0   =idum(4)
    nsteps=idum(5)
 print,idum
 readu,1,comega,calpha,reynolds,etadf,ybl,rzc,to,reyntp,rbot ; & etadf=0.001
 print,comega,calpha,reynolds,etadf,ybl,rzc,to ; & etadf=0.001
 nnodes=numnp/2
 etae=(2.6e-6)*(6.96e10)^2/comega/(1.e4)      ; envelope diffusivity (m^2/s)
 tdif=comega/(2.6e-6)/(3600.*24.*365.25)      ; diffusion time, in years
 u0=reynolds*(2.6e-6)*(6.96e10)/comega/1.e2   ; surface meridional flow (m/s)
 gam0=reyntp*(2.6e-6)*(6.96e10)/comega/1.e2 ; turbulent pumping speed (m/s)
 nt=nt0*nsteps

;---------- set some default values
if rbphi le 0.5 or rbphi gt 1.0 then rbphi=rzc
if latbphi le 0. or latbphi ge 90. then latbphi=45.

xx   =fltarr(nelx+1)
yy   =fltarr(nely+1)
b_r =fltarr(nt+1,nelx+1)
b_y =fltarr(nt+1,nelx+1)
b_yrcut =fltarr(nt+1,nely+1)
omg =fltarr(nelx+1,nely+1)
eta =fltarr(nelx+1,nely+1)
ur =fltarr(nelx+1,nely+1)
uth =fltarr(nelx+1,nely+1)
aquad=fltarr(nelx+1)
daqdm=fltarr(nelx+1)
bquad=fltarr(nelx+1,nely+1)
bphi=fltarr(nelx+1,nely+1)
time=findgen(nt+1)/(nt)*to ;*(1.*nt/nt0)
emag=fltarr(nt+1)

readu,1,xx                   ; latitudinal variable, mu=cos(theta)
readu,1,yy                   ; depth variable r/R
readu,1,omg
readu,1,ur
readu,1,uth
readu,1,eta
xxa=acos(xx)                 ; polar angle in radian...
xxan=-(2.*xxa/!PI-1.)*90.    ; ...here converted to latitude in degrees
print,yy(0:32)

;window,10,retain=2
;plot,yy,eta(nelx/2:nelx/2,*),xrange=[0.5,1.],psym=5
;oplot,yy,eta(0:0,*),linestyle=2
;oplot,yy,eta(nelx:nelx,*),linestyle=3
;oplot,[0.7,0.7],[0.,1.],linestyle=1
;oplot,[0.,1.],[0.5,0.5],linestyle=1
;
;window,11,retain=2
;plot,yy,omg(nelx/2:nelx/2,*),xrange=[0.5,1.],psym=5,yrange=[0.8,1.]
;oplot,[0.7,0.7],[0.,1.],linestyle=1
;---------- find nodal indices of desired cuts
ibp=0 & ibr=0 & jbp=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
for k=1,nelx-1 do if xxan(k+1) gt latbphi  and xxan(k) le latbphi then jbp=k
print, 'Bphi Butterfly diagram constructed at r/R=',yy(ibp)
print, 'Br   Butterfly diagram constructed at r/R=',yy(ibr)
print, 'Bphi Radial cut constructed at latitude=',xxan(jbp)

;-----------------------------------------------------------------------
; select contour interval and various plotting stuff
cin=[ 0,  10,     11,     12,     13,     14,     15,     16,     17,    $
           18,     20,         21,     22,     23,     24,     25,   $
           26,     27,    28,     29,     30, 31    ]+0
lev2=[-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,$
         0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1]
lev3=[-1.1,-1.0,-0.9,-0.8,-0.7,-0.5,-0.3,-0.2,-0.1,-0.01,-0.001,$
         0.001,0.01,0.1,0.2,0.3,0.5,0.6,0.7,0.8,0.9,1.0,1.1]
pos1 =[0.15,0.55,0.95,0.775]
pos2 =[0.15,0.3,0.95,0.525]
pos3 =[0.15,0.05,0.95,0.275]
a=2.*!PI*findgen(16)/15.
usersym,sin(a),cos(a),/fill
!X.STYLE=1 & !X.THICK=2
!Y.STYLE=1 & !Y.THICK=2

;---------- make up color table
device,decomposed=0
makecolortable

;---------- set up some labels
title1='!18B!D!9P!n!17(!18r/R=!17'+string(yy(ibp),format='(f4.2)')+'!17)'
title4='!18B!D!9P!n!17(!7h!17='+string(xxan(jbp),format='(f4.1)')+'!17)'
title2='!18B!Dr!n!17(!18r/R=!17'+string(yy(ibr),format='(f4.2)')+'!17)'
title3='!17'+title0
reycmlab='!18R!DM!n!17 =  '+string(reynolds,format='(f4.0)')
reytplab='!18R!DT!n!17 =  '+string(reyntp,format='(f4.0)')
wlab='!18w/R!17='+string(ybl,format='(f5.3)')
etalab='!7g!D!17c!n!17/!7g!D!170!n!17 ='+string(1./etadf,format='(f6.3)')
etaelab='!7g!D!17e!n!17 (m!U2!n s!U-1!n!17) ='+string(etae,format='(e8.1)')
colab='!18C!D!7X!n!17 ='+string(comega,format='(e8.1)')
calab='!18C!D!7a!n!17 ='+string(calpha,format='(f5.1)')
rbotlab='!18r!D!17b!n!17 ='+string(rbot,format='(f6.3)')
rzclab='!18r!D!17c!n!17 ='+string(rzc,format='(f6.3)')
tdiflab='!7s!17 (yr) ='+string(tdif,format='(e8.1)')
u0lab='!18u!D0!n!17 (m s!U-1!n!17) ='+string(u0,format='(e8.1)')
gam0lab='!7c!D0!n!17 (m s!U-1!n!17) ='+string(gam0,format='(e8.1)')
phaselab='!7D!9P!17/!7p!17 ='

;-----------------------------------------------------------------
; Read dynamo solution
istep=lonarr(1)
vdum=fltarr(5)
ifr=-1

for it=0,nt do begin    ; loop over time steps
   readu,1,istep
   readu,1,vdum & emag(it)=vdum(1)
   readu,1,bquad
   readu,1,bphi

;---------- Extract toroidal/poloidal fields at desired depths
    for j=0,nelx do begin
       b_y(it,j) =bphi(j,ibp)
       aquad(j)  =bquad(j,ibr)
    endfor

;---------- Extract toroidal field at selected latitude
    for k=0,nely do begin
       b_yrcut(it,k)=bphi(jbp,k)
    endfor

;---------- Compute vector potential angular derivative (to get B_r)
;       this computes d/dmu (A sin theta)    ( mu=cos(theta) )
    for j=1,nelx-1 do begin
       sinthp=sin(xxa(j+1))
       sinthm=sin(xxa(j-1))
       delmu=xx(j+1)-xx(j-1)
       daqdm(j)=(aquad(j+1)*sinthp-aquad(j-1)*sinthm)/delmu
    endfor
; patch for pole
    sinthp=sin(acos(xx(1)))
    sinthm=sin(acos(xx(0)))
    delmu=xx(1)-xx(0)
    daqdm(0)=(aquad(1)*sinthp-aquad(0)*sinthm)/delmu
; patch for equator
    sinthp=sin(acos(xx(nelx)))
    sinthm=sin(acos(xx(nelx-1)))
    delmu=xx(nelx)-xx(nelx-1)
    daqdm(nelx)=(aquad(nelx)*sinthp-aquad(nelx-1)*sinthm)/delmu

;---------- Now finally calculate B_r at surface
    for j=0,nelx do b_r(it,j)=-daqdm(j)

endfor       ; end of loop over time steps

;---------- Normalize fields to [-1.,1.], for plotting purposes
maxbr=max(b_r(nt/2:nt,*))
maxbyrcut=max(b_yrcut(nt/2:nt,*),irmax)
maxby=max(b_y(nt/2:nt,*),imax)
b_r    =b_r/maxbr 
b_y    =b_y/maxby
b_yrcut=b_yrcut/maxbyrcut
imax=jbp

;------------------------------------------------------------------------------
; Now a bit of Fourier analysis

fundamental=2.*!PI/to           ; fundamental frequency
nyquist=!PI/(to/nt)             ; Nyquist frequency
cutoff=fix(to/(1.*(to/nt)))     ; cutoff wavenumber
freq=findgen(cutoff)            ; frequencies, in units of fundamental frequency
spectrum=fft(b_y(*,imax))       ; fft of toroidal field timeseries
power=abs(spectrum)             ; power spectrum
peakpower=max(power(0:100),jj)   ; find frequency at peak power
peakfreq=freq(jj)*fundamental   ; peak frequency in spectrum
peakperiod=2.*!PI/peakfreq*tdif ; period (in yr) corresponding to peak freq.
periodlab='!18P !7[!17yr!7]!17 ='+string(peakperiod,format='(f6.2)')


;---------- set up more labels
maxbplab='!17max(!18B!D!9P!n!17) !7[!17T!7]!17 =' +string(maxby/1.e4,format='(e8.2)')
maxbrlab='!17max(!18B!Dr!n!17) !7[!17T!7]!17 ='   +string(maxbr/1.e4,format='(e8.2)')

;------------------------------------------------------------------------------
; now do the plots
; First window: time-latitude and radius-latitude diagram
window,1,retain=2,xs=700,ys=600
ti=to/2.     ; plot only second half of run

;---------- 1. plot color-filled contours for B-phi
    contour,pos=pos1,levels=lev2,/fill,$
       xcharsize=0.0000000001,charsize=1.25,c_color=cin,$
       xrange=[ti,to],yrange=[0.,90.],xticklen=0.04,$
       ytitle='!17Lat. (deg)',$
       title='!17 ',$
       b_y,time,xxan
    xyouts,charsize=1.5,to*1.007,90.*0.95,'!17NP'
    xyouts,charsize=1.5,to*1.007,0.0,'!17EQ'
    xyouts,ti-(to-ti)/10.,90.*0.2,orientation=90.,charsize=1.75,title1

    for ii=0,nt-1 do begin
       dum=max(abs(b_y(ii:ii,0:nelx)),jmax)
       oplot,psym=8,symsize=0.5,[time(ii),time(ii)],[xxan(jmax),xxan(jmax)],color=3
    endfor
    xyouts,ti+0.,90.*1.8,charsize=2.5,'!17'+solfile
;   xyouts,ti+0.,90.*1.6,charsize=2.5,'!17'+title0

    xyouts,ti-0.1*(to-ti),90.*1.6,charsize=1.75,calab
    xyouts,ti-0.1*(to-ti),90.*1.45,charsize=1.75,colab
    xyouts,ti-0.1*(to-ti),90.*1.3,charsize=1.75,reycmlab
    xyouts,ti-0.1*(to-ti),90.*1.15,charsize=1.75,reytplab

    xyouts,ti+0.15*(to-ti),90.*1.6,charsize=1.75,etalab
    xyouts,ti+0.15*(to-ti),90.*1.45,charsize=1.75,rzclab
    xyouts,ti+0.15*(to-ti),90.*1.3,charsize=1.75,rbotlab
    xyouts,ti+0.15*(to-ti),90.*1.15,charsize=1.75,wlab

    xyouts,ti+0.40*(to-ti),90.*1.6,charsize=1.75,color=5,etaelab
    xyouts,ti+0.40*(to-ti),90.*1.45,charsize=1.75,color=5,tdiflab
    xyouts,ti+0.40*(to-ti),90.*1.3,charsize=1.75,color=5,u0lab
    xyouts,ti+0.40*(to-ti),90.*1.15,charsize=1.75,color=5,gam0lab

    xyouts,ti+0.75*(to-ti),90.*1.6,charsize=1.75,color=6,maxbplab
    xyouts,ti+0.75*(to-ti),90.*1.45,charsize=1.75,color=6,maxbrlab
    xyouts,ti+0.75*(to-ti),90.*1.3,charsize=1.75,color=6,periodlab
;   xyouts,ti+0.75*(to-ti),90.*1.15,charsize=1.75,color=6,phaselab

    contour,/noerase,pos=pos2,levels=lev2,/fill,$
       xcharsize=0.0000000001,charsize=1.25,c_color=cin,$
       xrange=[ti,to],yrange=[0.5,1.],xticklen=0.04,$
       ytitle='!18r/R!17',$
       title='!17 ',$
       b_yrcut,time,yy
    oplot,linestyle=5,[ti,to],[rzc,rzc],thick=2
    xyouts,charsize=1.5,to*1.007,1.0*0.98,'!17TOP'
    xyouts,charsize=1.5,to*1.007,rzc*0.98,'!17ZC'
    xyouts,ti-(to-ti)/10.,0.6+0.2*0.4,orientation=90.,charsize=1.75,title4
    
    colorscale,pos2,lev2,cin,-0.18,-0.13,-0.5,1.5,2

    contour,/noerase,pos=pos3,levels=lev3,/fill,$
       c_linestyle=clst3,charsize=1.25,c_color=cin,$
       xrange=[ti,to],yrange=[0.,90],xticklen=0.04,$
       ytitle='!6Lat. (deg)',$
       xtitle='!18t/!7s!17',$
       b_r,time,xxan
    xyouts,charsize=1.5,to*1.007,90.*0.95,'!17NP'
    xyouts,charsize=1.5,to*1.007,0.0,'!17EQ'
    xyouts,ti-(to-ti)/10.,90.*0.2,orientation=90.,charsize=1.75,title2

;------------------------------------------------------------------------------
; Second window: time series of toroidal field, and power spectrum


window,2,retain=2,xs=500,ys=700
nf=100.
yr=[1.e-4,1.]
plot_io,pos=[0.175,0.4,0.95,0.95],freq(0:nf),power(0:nf),psym=-8,$
   yrange=yr,charsize=1.5,xtitle='!7x!D!18k!n!17/!7x!D!170!n!17',$
   ytitle='!17Power',title='!17'+solfile
oplot,linestyle=1,[freq(jj),freq(jj)],yr
oplot,linestyle=1,[2.*freq(jj),2.*freq(jj)],[yr(0),yr(0)*10.]
oplot,linestyle=1,[3.*freq(jj),3.*freq(jj)],[yr(0),yr(0)*10.]
oplot,linestyle=1,[4.*freq(jj),4.*freq(jj)],[yr(0),yr(0)*10.]
oplot,linestyle=1,[5.*freq(jj),5.*freq(jj)],[yr(0),yr(0)*10.]

;----- some labels and numerical info
delx=0.25*(freq(nf)-freq(0))
freqlab='!17Peak frequency (!7s!U!17-1!n!17) : '+string(peakfreq,format='(e10.3)')
periodlab='!17Peak period (yr) : '+string(peakperiod,format='(e10.3)')
xyouts,charsize=1.5,delx,0.5,freqlab
xyouts,charsize=1.5,delx,0.3,periodlab

;----- the original time series from which FFT was computed
ytit=$
'!18B!D!9P!n!17(!18r/R=!17'+string(yy(ibp),format='(f4.2)')+$
'!17,!7h!17='+string(xxan(imax),format='(f4.1)')+'!17)'
plot,xrange=[0.,to],pos=[0.175,0.1,0.95,0.3],/noerase,time,b_y(*,imax),$
       xtitle='!18t/!7s!17',thick=2,yrange=[-1.,1.],charsize=1.5,$
       ytitle=ytit

return
end