FUNCTION read_input, io, delim, param
   line = ' '
   WHILE ((NOT STRMATCH(line,'*'+param+'*')) AND (NOT EOF(io))) DO READF,io,line
   IF (EOF(io)) THEN BEGIN
      PRINT, param, ' not found in run.out'
      POINT_LUN, io, 0
      RETURN, '0'
   ENDIF
   indices = STRSPLIT(line, delim)
   string = STRMID(line,indices[1],20)
   RETURN, string
END

PRO gyro_dens_rtheta_out, sim_dir, THETA_MAX=theta_max

; Inputs: 	sim_dir - GYRO simulation directory containing the files
;		   t.out, profile_vugyro.out, run.out, profile_sim.out, and moment_n.out.
;		   If geometry_arrays.out exists, then input.dat is also necessary.
;		   If run from simdir, set to ''. Otherwise must end with '/'.
;
;  		theta_max - Extent of poloidal domain (rad)
;		   (theta_min=-theta_max)
;		   If left off, extent will be one wavelength of longest mode:
;		   L_pol = 2*pi*rnorm/n_min/q.
; Output:  	densities = densities in units of 10^19 m^-3 in region
;  		   symmetrical about outer midplane,
;  		   written to 'densities_rtheta_out.sav' in sim_dir,
;  		   array of length (n_s,n_th,n_r,n_t)
;  		   where n_s = number of species
;  			densities(0,...) = main ion
;  			densities(1,...) = imp. ion 1
;  			densities(2,...) = imp. ion 2
;  			densities(n_s-1,...) = electron
;  		   n_th = number of theta grid points
;  			= FIX(n_r*theta_max*rnorm/L_r),
;  			  where L_r is the radial width of the domain.
;  		   n_r = number of radial grid points (= RADIAL_GRID)
;  		   n_t = number of time points
;               rmin = midplane minor radius (m), of length n_r
;  		theta = poloidal grid (rad) of length n_th
;  		   (zero is at outer midplane)
;  		t = time grid (a/cs) of length n_t
;  		kappa = plasma elongation, of length n_r
;  		delta_0 = average plasma triangularity, of length n_r
;  		        = (delta_up+delta_low)/2
;  		delta_1 = differential plasma triangularity, of length n_r
;  		        = (delta_up-delta_low)/2
;		R0 = major radius (m) of center of flux surfaces defined by rmin,
;		     of length n_r
;		z0 = elevation (m) of center of flux surfaces defined by rmin,
;		     of length n_r
;		a = rmin at separatrix (m)
;		csda = cs/a
;		omega_eb_0 = ExB Doppler frequency for n = 1
;		             at the norm point of the box
;
; Notes:  The GYRO input parameter THETA_PLOT must be even.
;  	  Otherwise, contact me and I'll modify the routine.
;
;  	  The density distributions can be constructed in (R,z) space
;	  using the following formulae:
;
;	  R(rmin,theta) = R0(rmin) + rmin*cos(theta+arcsin(delta(rmin))*sin(theta))
;			where delta(rmin) = delta_0(rmin)+delta_1(rmin) for theta > 0
;					  = delta_0(rmin)-delta_1(rmin) for theta < 0
;	  z(rmin,theta) = z0(rmin) + kappa(rmin)*rmin*sin(theta)
;
; Revisions:
;  	8/21/05:  Made compatible with global runs.
;  	 3/1/06:  Ditto for GYRO v. 5.0.x
;       3/23/06:  Added outputs rho, R_maj
;       2/17/07:  Removed outputs rho, R_maj
;		  Added outputs kappa, delta_0, delta_1, R0, z0,
;		  a, csda, omega_eb_0
;       5/17/07:  Added capability to read nu(r,theta) to replace q*theta

pi = !PI

PRINT, 'Reading parameters from run.out ...'
OPENR, 1, sim_dir+'run.out', ERROR = err
   IF (err NE 0) THEN BEGIN
      PRINT, !ERROR_STATE.MSG
      RETURN
   ENDIF
   rp_method = read_input(1, ':', 'profile_method')
      IF (STRMATCH(rp_method, '*EXP*')) THEN rp_method=3 ELSE rp_method=5
   n_n = FIX(read_input(1, ':', 'n_n'))
   n_r = FIX(read_input(1, ':', 'n_x'))
   n_kin = FIX(read_input(1, ':', 'n_kinetic'))
   radius = FLOAT(read_input(1, ':', '.RADIUS'))
   Roa = FLOAT(read_input(1, ':', 'ASPECT_RATIO'))
   q0 = FLOAT(read_input(1, ':', 'SAFETY_FACTOR'))
   s_hat = FLOAT(read_input(1, ':', 'SHEAR'))
   rho_star = FLOAT(read_input(1, ':', 'RHO_STAR'))
   nione = FLOAT(read_input(1, ':', 'NI_OVER_NE'))
   IF (n_kin GE 3) THEN nz1one = FLOAT(read_input(1, ':', 'NI_OVER_NE_2'))
   IF (n_kin GE 4) THEN nz2one = FLOAT(read_input(1, ':', 'NI_OVER_NE_3'))
   betae_unit = FLOAT(read_input(1, ':', 'BETAE_UNIT'))
   mu_e = FLOAT(read_input(1, ':', 'MU_ELECTRON'))
;  For versions earlier than 5.2.x:
   IF (Roa LT 0.0001) THEN BEGIN
      Ror = FLOAT(read_input(1, ':', 'R_major'))
      Roa = Ror*radius
   ENDIF
;
   omega_eb_0 = FLOAT(read_input(1, ':', 'omega_eb_0'))
   IF (omega_eb_0 EQ 0) THEN BEGIN
      PRINT, '  Assuming omega_eb_0 = 0'
      omega_eb_0 = 0.
   ENDIF
   a = FLOAT(read_input(1, ':', 'a_unit_norm (m)'))
      rnorm = radius*a  ;  radial center of box (m)
   csda = FLOAT(read_input(1, ':', 'csda_norm'))
   ne0 = FLOAT(read_input(1, ':', 'den_norm'))
   line = ' '
   WHILE (NOT STRMATCH(line, '*WAVENUMBERS*')) DO READF, 1, line
   n0 = FIX(read_input(1, '=', 'n'))
   n1 = FIX(read_input(1, '=', 'n'))
   n_dn = n1 - n0
CLOSE, 1
PRINT, '  Done'

PRINT, 'Reading output times from t.out ...'
OPENR, 2, sim_dir+'t.out', ERROR = err
   IF (err NE 0) THEN BEGIN
      PRINT, !ERROR_STATE.MSG
      RETURN
   ENDIF
   i = 0
   WHILE (NOT EOF(2)) DO BEGIN
       READF, 2, dummy
       i = i+1
   ENDWHILE
CLOSE, 2
n_t = i
tt = FLTARR(2,n_t)
OPENR, 2, sim_dir+'t.out'
   READF, 2, tt
CLOSE, 2
PRINT, '  Done'
t = FLTARR(n_t)
t[*] = tt[1,*]

;;--------------------------Only for testing----------------------
;n_t = 2
;;----------------------------------------------------------------

rmin = FLTARR(n_r)
q = FLTARR(n_r)
kappa = FLTARR(n_r)
delta_0 = FLTARR(n_r)
delta_1 = FLTARR(n_r)
R0 = FLTARR(n_r)
z0 = FLTARR(n_r)

PRINT, 'Reading data from profile_vugyro.out ...'
OPENR, 4, sim_dir+'profile_vugyro.out', ERROR = err
   IF (err NE 0) THEN BEGIN
      PRINT, !ERROR_STATE.MSG
      RETURN
   ENDIF
   ;  Read value of THETA_PLOT:
   n_theta = 1
   SKIP_LUN, 4, 5, /LINES
   READF, 4, n_theta
   PRINT, 'THETA_PLOT = ', n_theta
   IF (n_theta EQ 1) THEN BEGIN
      PRINT, '  You are requesting densities at the outer midplane when'
      PRINT, '  coefficients are provided only at the inner midplane.'
      PRINT, '  Exiting ...
      RETURN
   ENDIF
   ;  Look for first possible value of rho:
   dum = '0'
   WHILE (STRLEN(STRCOMPRESS(dum)) LT 6) DO READF, 4, dum
   rmin[0] = FLOAT(dum)
   FOR i=1,n_r-1 DO BEGIN
      READF, 4, dum
      rmin[i] = dum
   ENDFOR
   rmin = a*rmin
;CASE rp_method OF
;  3: BEGIN
   dens_norm = FLTARR(n_kin,n_r)
   READF, 4, q
   SKIP_LUN, 4, 2*n_r, /LINES
   SKIP_LUN, 4, 3*n_kin*n_r, /LINES
   READF, 4, dens_norm
   SKIP_LUN, 4, 2*n_r, /LINES
   READF, 4, delta_0
   READF, 4, delta_1
   READF, 4, kappa
CLOSE, 4
;      END
;  5: BEGIN
;	q = q0*(s_hat*(rho/radius-1.) + 1.)
;	A_ref = mu_e^2/1836.  ;  atomic mass of main ion
;	rho_s = rho_star*a  ;  ion gyroradius at electron temperature (m)
;	ne0 = betae_unit/rho_s^2*A_ref/387.4  ;     average n_e (10^19 m^-3)
;	r = a*rho
;	CLOSE, 4
;      END
;   ELSE: BEGIN
;	PRINT, 'Invalid RADIAL_PROFILE_METHOD'
;	RETURN
;      END
;ENDCASE
PRINT, '  Done'

;  Obtain R0 and z0:

PRINT, 'Reading data from profile_sim.out ...'
OPENR, 5, sim_dir+'profile_sim.out', ERROR = err
   IF (err NE 0) THEN BEGIN
      PRINT, !ERROR_STATE.MSG
      RETURN
   ENDIF

;   SKIP_LUN, 5, (n_r+4)*(n_kin+1)+2, /LINES
;   (Removed by RVB, 5/17/07)  Replaced by:
   line = ' '
   WHILE ((NOT STRMATCH(line,'*aspect*')) AND (NOT EOF(5))) DO READF,5,line
   IF (EOF(5)) THEN BEGIN
      PRINT,  'aspect_s not found in profile_sim.out'
      CLOSE, 5
      RETURN
   ENDIF
   SKIP_LUN, 5, 1, /LINES
   FOR i=0,n_r-1 DO BEGIN
      READF, 5, FORMAT='(12x,G11.4)', aspect_s
      R0[i] = aspect_s*rmin[i]
   ENDFOR

;   SKIP_LUN, 5, (n_r+4)*3+2, /LINES
;   (Removed by RVB, 5/17/07)  Replaced by:
   line = ' '
   WHILE ((NOT STRMATCH(line,'*z_mag*')) AND (NOT EOF(5))) DO READF,5,line
   IF (EOF(5)) THEN BEGIN
      PRINT,  'z_mag_axis not found in profile_sim.out'
      PRINT,  ' setting to zero ...'
      z0 = 0.*z0
   ENDIF ELSE BEGIN
      SKIP_LUN, 5, 1, /LINES
      FOR i=0,n_r-1 DO BEGIN
	READF, 5, FORMAT='(48x,G11.4)', z_mag
	z0[i] = z_mag
      ENDFOR
   ENDELSE
CLOSE, 5
PRINT, '  Done'

;; Dimensions of box:
L_r = rmin[n_r-1] - rmin[0]

;  All new by RVB, 5/17/07
;------------------------------------------------------------------
IF ~ KEYWORD_SET(theta_max) THEN theta_max=pi/n_dn/q0
L_pol = theta_max*rnorm

n_th = FIX(n_r*L_pol/L_r)
nu = FLTARR(n_r,n_th)
theta_min = -theta_max
theta = theta_min+findgen(n_th)*(theta_max-theta_min)/(n_th-1)

OPENR, 6, sim_dir+'geometry_arrays.out', ERROR = err
  IF (err EQ 0) THEN BEGIN
;    Read THETA_MULT:
     OPENR, 7, sim_dir+'input.dat', ERROR = err
     IF (err NE 0) THEN BEGIN
	PRINT, !ERROR_STATE.MSG
	CLOSE, 6
	RETURN
     ENDIF
     line = ' '
     WHILE ((NOT STRMATCH(line,'*THETA_MULT*')) AND (NOT EOF(7))) DO READF,7,line
     IF (EOF(7)) THEN BEGIN
        PRINT,  'THETA_MULT not found in input.dat'
	CLOSE, 7
        RETURN
     ENDIF
     CLOSE, 7
     indices = STRSPLIT(line, ' ',/EXTRACT)
     theta_mult = FIX(indices[0])
     PRINT, 'THETA_MULT = ', theta_mult
     PRINT,'Reading nu(r,theta)  from geometry_arrays.out ...'
     n_th_nu = n_theta*theta_mult
     nu_in = FLTARR(n_r,n_th_nu)
     FOR i = 0, n_r-1 DO BEGIN
	FOR j = 0, n_th_nu-1 DO BEGIN
	   READF, 6, temp
	   nu_in[i,j] = temp
	ENDFOR
     ENDFOR
     theta_in = -pi + 2*pi/n_th_nu*FINDGEN(n_th_nu)
     FOR i=0,n_r-1 DO BEGIN
 	glop = REFORM(nu_in[i, *])
	nu[i,*] = INTERPOL(glop, theta_in, theta)
     ENDFOR
     PRINT, '  Done'
  ENDIF ELSE BEGIN
     FOR i=0,n_r-1 DO nu[i,*] = q[i]*theta
  ENDELSE
CLOSE, 6
;------------------------------------------------------------------

;; toroidal mode number
n = n0+indgen(n_n)*n_dn

;;  For n_theta even:

theta0 = 0.
i0 = FIX(n_theta/2)
thetabot = -pi/i0
ibot = i0 - 1
thetatop = pi/i0
IF (n_theta ne 2) THEN $
  itop = i0 + 1 ELSE $
  itop = i0 - 1

c_n = FLTARR(2,n_theta,n_r,n_kin,n_n,n_t)
c = complexarr(n_theta,n_r,n_kin,n_n)
cc = complexarr(n_kin,n_r,n_th)

;;-----------------------------
PRINT, 'Reading data from moment_n.out ...'
OPENR, 5, sim_dir+'moment_n.out'
READF, 5, c_n
CLOSE, 5
PRINT, '  Done'
;;-----------------------------

PRINT, 'Calculating density array ...'

density = FLTARR(n_kin,n_r,n_th,n_t)

;; step through time
FOR it=0,n_t-1 DO BEGIN

  PRINT, '  time = ', t(it)
  ;; complex, slowly-varying density coefficients at selected time c[theta,r,s,n]:
  c[*,*,*,*] = complex(c_n[0,*,*,*,*,it],c_n[1,*,*,*,*,it])

  ;; n=0
  FOR is = 0,n_kin-1 DO BEGIN
    FOR j=0,n_th-1 DO BEGIN
     IF (theta[j] LT theta0) THEN BEGIN
        cc[is,*,j] = c[ibot,*,is,0]*(theta0-theta[j])/(theta0-thetabot)+$
                  c[i0,*,is,0]*(theta[j]-thetabot)/(theta0-thetabot)
     ENDIF ELSE BEGIN
        cc[is,*,j] = c[itop,*,is,0]*(theta[j]-theta0)/(thetatop-theta0)+$
                  c[i0,*,is,0]*(thetatop-theta[j])/(thetatop-theta0)
     ENDELSE
    ENDFOR
    density[is,*,*,it] = FLOAT(cc[is,*,*])
  ENDFOR


  ;; n>0
  FOR is = 0,n_kin-1 DO BEGIN
    FOR in=1,n_n-1 DO BEGIN
      FOR j=0,n_th-1 DO BEGIN
        IF (theta[j] LT theta0) THEN BEGIN
          cc[is,*,j] = c[ibot,*,is,in]*(theta0-theta[j])/(theta0-thetabot)+$
                  c[i0,*,is,in]*(theta[j]-thetabot)/(theta0-thetabot)
        ENDIF ELSE BEGIN
          cc[is,*,j] = c[itop,*,is,in]*(theta[j]-theta0)/(thetatop-theta0)+$
                  c[i0,*,is,in]*(thetatop-theta[j])/(thetatop-theta0)
        ENDELSE

        density[is,*,j,it] = density[is,*,j,it] + 2*FLOAT(cc[is,*,j]* $
                               exp(-complex(0,1)*n[in]*nu[*,j])* $
                               exp(-complex(0,1)*n[in]*omega_eb_0*t[it]))
      ENDFOR
    ENDFOR
    IF (rp_method EQ 5) THEN BEGIN
	IF (is EQ 0) THEN dens0 = nione
	IF (is EQ 1) AND (n_kin GT 2) THEN dens0 = nz1one
	IF (is EQ 2) AND (n_kin GT 3) THEN dens0 = nz2one
	IF (is EQ n_kin-1) THEN dens0 = 1.
	density[is,*,*,it] = ne0*(density[is,*,*,it] + dens0)
    ENDIF ELSE BEGIN
	FOR ix=0,n_r-1 DO BEGIN
	    density[is,ix,*,it] = ne0*(density[is,ix,*,it] + dens_norm[is,ix])
	ENDFOR
    ENDELSE
  ENDFOR

ENDFOR

PRINT, '  Done'
PRINT, 'Saving to "densities_rtheta_out.sav"'

save, filename=sim_dir+'densities_rtheta_out.sav', density, rmin, theta, t, $
	kappa, delta_0, delta_1, R0, z0, a, csda, omega_eb_0

END