; FILE: helio14_anim.idl  [IDL source code]
; production mpeg selon la méthode-Patrice
;===========================================================
; NONLINEAR AXISYMMETRIC BABCOCK-LEIGHTON DYNAMO
;-----------------------------------------------------------
; APAS7500 Fall 1997  Problem II.X.Y
; see Section X.Y of notes, and Appendix X
;-----------------------------------------------------------
; before running this first compile the following IDL procedures:
;    colorscale_pro.idl
;    mkcolort_pro.idl
;-----------------------------------------------------------
;-----------------------------------------------------------
; INPUT: this code requires an input file, the name of which
; being set immediately below
;===========================================================
 imakemovie=1
 solfile='saasfe/test_corinne.i3e' 
;solfile='saasfe/BL_ivp.ref2.i3e' 
 solfile='PLL-5.i3e' 
 titlelab='!17BL C!D!7a!n!17=+5  !7Dg=1/300 !17Rm=840' 
;imakemovie=0

 parity=-1.
;parity=0.
;device,decomposed=0
;makecolortable
;tvlct,rr,gg,bb,/get
 a=findgen(16)*(!PI*2./16.)
 usersym,cos(a),sin(a),/fill
 !ORDER=1
 imovie=1
;tvlct,red,green,blue

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

 title0=string(' ',format='(a80)') & ctime=1.e0
 readu,1,title0  & titlelab=title0
 print,'Using Input file: ',solfile
 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 ; & etadf=0.001
 print,comega,calpha,reynolds,etadf,ybl,rzc,to ; & etadf=0.001
 nt0=nt0*nsteps
 nnodes=numnp/2
;nskip=1500 ; number of time steps skipped over; for sim19hr
 nskip= 500 ; number of time steps skipped over; for sim6 et sim9
;nskip=0 ; number of time steps skipped over
 nt=500    ; number fo time steps in animation
;if nt0 lt nskip+nt then nt=max(nt0-nskip,1)
nt=100

xx   =fltarr(nelx+1)
yy   =fltarr(nely+1)
b_r =fltarr(nt+1,nelx+1)
b_y =fltarr(nt+1,nelx+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)
daqdt=fltarr(nelx+1)
bquad=fltarr(nelx+1,nely+1)
b_y=fltarr(nelx+1,nely+1)
to=to*nt/nt0
time=findgen(nt+1)/(nt)*to ;*(1.*nt/nt0)
emag=fltarr(nt+1)

readu,1,xx
readu,1,yy
readu,1,omg
readu,1,ur
readu,1,uth
readu,1,eta
xxa=acos(xx)

yc   =fltarr(nelx+1,nely+1)
xc   =fltarr(nelx+1,nely+1)

;--------- construct 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
;---------- find boundaries of layer and top
iblay=0
ilay=0
iint=0
for k=1,nely-1 do if yy(k+1) gt rzc  and yy(k) le rzc  then iint=k
for k=1,nely-1 do if yy(k+1) gt (rzc-ybl) and yy(k) le (rzc-ybl) then iblay=k
for k=1,nely-1 do if yy(k+1) gt 1.0 and yy(k) le 1.0 then itop=k
print, iblay,yy(iblay)
print, iint,yy(iint)
print, itop,yy(itop)

;-----------------------------------------------------------------------
; select contour interval and color levels

;   1. choose contour interval and color levels
 cin=[ 0,  10,     11,     12,     13,     14,     15,     16,     17,    $
 	   18,     19,     20,     21,     22,     23,     24,     25,   $
 	   26,     27,     28,     29,     30    ]+1
levb_y=[-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]
levquad=[-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.3,-0.1,-0.03,-0.01,$
         0.01,0.03,0.1,0.3,0.5,0.6,0.7,0.8,0.9,1.0,1.1]
;levquad=levb_y
;-----------------------------------------------------------------------------
; prepare output file for movie

 if imovie eq 1 then begin ; workstation animation
    nnx=350
    nny=350
    pos0=[0.1,0.1,0.975,0.975]
    mimage=bytarr(nnx,nny,nt+1)
    itr=0
    xr=[0.,1.1]
    yr=[0.,1.1]
    xlab=0.3
    ylab=-0.1
    xlab2=0.
    ylab2=-0.1
    thk=2
    chsiz=1.5
    ncyc=1.
    ccol=   [7,7,7,7,7,7,7,7,7,7,7,3,3,3,3,3,3,3,3,3,3,3]
 endif
;   nnx=275
;   nny=275
 window,0,retain=2,xs=nnx,ys=nny
 xinteranimate,set=[nnx,nny,nt],/showload,mpeg_bitrate=104857200.0,/mpeg_open,mpeg_filename='PLL-5.mpg'

;device,decomposed=0
makecolortable,red,green,blue
tvlct,red,green,blue

relab='!6Re(s)'
imlab='!6Im(s)'
colab='!6C!7X!6'
calab='!6C!7a!6'
;relab2='!6= '+string(reomega,format='(f8.3)')
;imlab2='!6= '+string(imomega,format='(f8.3)')
colab2='!6=1.E5'
calab2='!6='+string(calpha,format='(f4.1)')

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

 istep=lonarr(1)
 vdum=fltarr(5)
 ifr=-1

;---------- read results for eigenvector and eigenvalue
;           (skip over nskip solutions)
 aqmax=0.
 bymax=0.
 for iskip=0,nskip-1 do begin

    readu,1,istep
    readu,1,vdum
    print,istep
    readu,1,bquad
    readu,1,b_y

    if iskip gt nskip/2 then begin
       aqtmax=max(abs(bquad)) & if aqtmax ge aqmax then aqmax=aqtmax
       bytmax=max(abs(b_y))   & if bytmax ge bymax then bymax=bytmax
    endif
 endfor

; read second half of sim
 for it=0,nt-1 do begin
;for it=0,10 do begin

    readu,1,istep
    readu,1,vdum
       time(it)=vdum(0) & emag(it)=vdum(2)
    print,istep,time(it)
    readu,1,bquad
    readu,1,b_y

    bquad=bquad/aqmax
    b_y=b_y/bymax

; plot results
;------------------------------------------------------------------------------
    !X.STYLE=5
    !Y.STYLE=5
;---------- 1. plot color-filled contours for B-phi
;   if ifr eq 0 then erase
    contour,pos=pos0,levels=levb_y,$
       /fill,c_color=cin,$
       xrange=xr,yrange=yr,$
       b_y(itr:nelx,0:itop),xc(itr:nelx,0:itop), yc(itr:nelx,0:itop)
   oplot,/noclip,xc(itr:nelx,itop:itop), yc(itr:nelx,itop:itop),thick=thk
   oplot,/noclip,xc(itr:nelx,iint:iint), yc(itr:nelx,iint:iint),thick=2,$
      linestyle=2
   oplot,/noclip,[0.,0.],[0.,1.],thick=thk
   oplot,/noclip,[0.,1.],[0.,0.],thick=thk
;---------- 2. overplot contours for B-pol
   contour,/follow,/noclip,/noerase,pos=pos0,levels=levquad,$
      xrange=xr,yrange=yr,$
      c_labels=[0,0,0,0,0,0,0,0,0],c_color=ccol,c_thick=thk,$ 
      bquad(itr:nelx,0:itop+8),xc(itr:nelx,0:itop+8), yc(itr:nelx,0:itop+8)
    steplab='!8k!6='+string(istep,format='(i4)')
    timelab='!8t/!7s!17='+string(time(it),format='(f6.4)')
    xyouts,xlab2,ylab ,charsize=chsiz,timelab,color=5
;   xyouts,charsize=1.2,-0.02,1.05,titlelab
    xyouts,charsize=1.5,0.05,0.05,'!17'+solfile
;   endif

;   5. color scale
    colorscale,pos0,levb_y,cin,-0.1,-0.05,0.,1.0,thk

;   6. store image
;  ifr=ifr+1
;  if imovie eq 1 then begin
;     mimage(*,*,ifr:ifr)=tvrd(0,0,nnx,nny)
;  endif

; -- store animation as mpeg file (fenetre de taille nnx x nny )
      xinteranimate, frame=it, window=0
;     erase

endfor ; fin boucle temporelle, compteur "it"

      xinteranimate,30,/keep_pixmaps,/mpeg_close

end