;note the cubes are written in the way that the second index is the vertical direction
;e.g. vyy(x,y,z) - x,z are horizontal directions, y is vertical direction, vyy is velocity in vertical direction

;the code stores conservative variables. Instead of velocity there is momentum. To get velocity the momentum has to be
;divided by the density

;also, the code stores the magnetic field in its internal units. The units are chosen in that way that the magnetic
;pressure is B^2/2. To convert the magnetic field components to cgs they have to be multiplied by sqrt(4*pi).



;path to models
ptk='./'

;model index
nmod='088200'


;model size horizontal directions
npt=288

;model size vertical direction
npty=100

close,1

;coefficient to convert internal code units into cgs units
coef=sqrt(4*3.14159265)


;get temperature and gas pressure (tem and pre)
openr,1,ptk+'eos.'+nmod
a=assoc(1,fltarr(npt,100,npt))
tem=a(0)
pre=a(1)
close,1
print,tem(0,0,0)


;get density rho
openr,1,ptk+'result_0.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
rho=a(0)
close,1

;get horizontal x component of momentum vxx
openr,1,ptk+'result_1.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
vxx=a(0)
close,1


;get vertical y component of momentum vyy
openr,1,ptk+'result_2.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
vyy=a(0)
close,1


;get horizontal z component of momentum vzz
openr,1,ptk+'result_3.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
vzz=a(0)
close,1


;get horizontal x component of magnetic field vector bxx
openr,1,ptk+'result_5.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
bxx=a(0)
close,1


;get vertical y component of magnetic field vector byy
openr,1,ptk+'result_6.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
byy=a(0)
close,1

;get horizontal z component of magnetic field vector bzz
openr,1,ptk+'result_7.'+nmod
a=assoc(1,fltarr(npt,npty,npt))
bzz=a(0)
close,1


;convert momenta to velocities
vxx=vxx/rho
vyy=vyy/rho
vzz=vzz/rho

;convert magnetic field components to cgs units
bxx=bxx*coef
byy=byy*coef
bzz=bzz*coef

tvframe,byy(*,60,*),/ba,/sa,/as



end