pro show_collapse
  for i=0,18 do begin
     snr=string(i,form='(i03)')
     readnew,"snap_"+snr,x,"POS"
     plot,x[0,*],x[1,*],psym=3,xtit='x',ytit='y',chars=2,charth=2,xth=2,yth=2,xr=[-1,1],yr=[-1,1],xst=1,yst=1
     save_screen,"frame_"+snr
  end
end

pro show_energies
  readcol,'energy.txt',comment='#',delimiter=' ',$
          format='F,F,F,F,*',$
          t,Eint,Epot,Ekin

  M = 1.5815
  R = 1.0
  G = 1.0
  eps_star = G * M / R
  rho_star = 3 * M / (4 * !Pi * R^3)
  p_star = eps_star * rho_star
  t_star = sqrt(R^3 / (G * M))
  print,' =>  M=',m,', rho_*=',rho_star,', p_* =',p_star,', t_* =',t_star,' eps_* =',eps_star

  plot,[1],[1],xr=[0,3],yr=[-3,2],xst=1,yst=1,xtit="t/t*",ytit="E/e*",chars=2,charth=2,xth=2,yth=2
  oplot,t/t_star,Eint/eps_star/M
  oplot,t/t_star,Epot/eps_star/M
  oplot,t/t_star,Ekin/eps_star/M
  oplot,t/t_star,(Eint+Epot+Ekin)/eps_star/M,l=2

end

pro show_profiles
   M = 1.5815
   R = 1.0
   G = 1.0
   eps_star = G * M / R
   rho_star = 3 * M / (4 * !Pi * R^3)
   p_star = eps_star * rho_star
   t_star = sqrt(R^3 / (G * M))
   print,' =>  M=',m,', rho_*=',rho_star,', p_* =',p_star,', t_* =',t_star,' eps_* =',eps_star
   gam1 = (5./3.-1)

   FOR k=1,24 DO BEGIN
     snr=string(k,form='(i02)')
  
     name_ppm='PPM/PPM_'+snr

     name='snap_0'+snr
     readnew,name,h,'HEAD'
     readnew,name,x,'POS'
     readnew,name,rho,'RHO'
     readnew,name,u,'U'
     readnew,name,m,'MASS'
     readnew,name,hsml,'HSML'

     tit='t/t*='+string(h.time/t_star,form='(f4.2)')
     p = gam1 * u * rho 

     plot_oo,[1],[1],/nodata,tit=tit,xr=[0.0002,2],yr=[0.04,4000],xst=1,yst=1,xtit='r/R',ytit='rho/rho* , p/p*/10',charsize=2,charthick=2

     rr = sqrt((x(0,*))^2+(x(1,*))^2+(x(2,*))^2)
     oplot,rr,rho/rho_star,psym=3
     oplot,rr,p/p_star/10,psym=3

     oplot,min(hsml)*[1,1],[0.001,1e6],col=64
   
     readcol,name_ppm,comment='#',delimiter=' ',format='F,F,F,F',r,dens,presh,dummy
     oplot,r,dens,col=128,thick=2
     oplot,r,presh/10,col=128,thick=2

     save_screen,"profile_0"+snr
  END

end