;+
;  $Id: load_dens.pro,v 1.2 2008/09/23 14:12:43 thernis Exp $
; PURPOSE:
;  Load a PFSS density file generated by Y-M Wang's program
; CATEGORY:
;  raytracing, io
; INPUTS:
;  file: filename of the density cube
; OUTPUTS:
;  dens: the density cube
;  rco: r axis parameters
;  phico: phi axis parameters
;  thetaco: theta axis parameters
;  r: r axis array
;  phi: phi axis array
;  theta: theta axis array
;
;-
pro load_dens,dens,rco,phico,thetaco,r,phi,theta,file=file,savefile=savefile
if n_elements(file) eq 0 then begin
    print,'Using default density file'
    rtinitenv
    file=getenv('RT_PATH')+'/data/scraytrace/pfssdt.fits'
endif
cube=mrdfits(file,1,/silent)
dh=(size(cube.dens))[1]
dl=(size(cube.dens))[2]
dt=(size(cube.dens))[3]

r=findgen(dh)*0.01+1.
rco=[1.,0.01,dh]
phi=findgen(dl)*(2*!pi)/float(dl)+(2*!pi)/180.
phico=[(2*!pi)/180.,(2*!pi)/float(dl),dl]
theta=-!pi/2.+!pi/90+!pi/dt*findgen(dt)
thetaco=[-!pi/2.+!pi/90,!pi/dt,dt]
dens=cube.dens

end