openr,1,'katrina_files.txt'
	
;	openw,3,'wsr_88d_output'
		
	file = ''					;define one string line

;	using 9 reflectivity sweeps
	car_r = fltarr(481,481,9)
	car_z = fltarr(481,481,9)
	
	cnt = -1

	FOR k = 0,20-1 DO BEGIN
	
	readf,1,file

	strings = strsplit(file,',', /EXTRACT)		;break up the string

	file = strings[0]
	sd = strings[1]
	sh = strings[2]
	sm = strings[3]
	ss = strings[4]


        id = ncdf_open(file)
        id1 = ncdf_varid(id,"Range_to_First_Cell")
        ncdf_varget, id, id1, range2fc
        id1 = ncdf_varid(id,"Cell_Spacing")
        ncdf_varget, id, id1, cellsp
        id1 = ncdf_varid(id,"Latitude")
        ncdf_varget, id, id1, lat
        id1 = ncdf_varid(id,"Longitude")
        ncdf_varget, id, id1, lon
        id1 = ncdf_varid(id,"Altitude")
        ncdf_varget, id, id1, alt
        id1 = ncdf_varid(id,"Azimuth")
        ncdf_varget, id, id1, azi
        id1 = ncdf_varid(id,"Elevation")
        ncdf_varget, id, id1, ele
        id1 = ncdf_varid(id,"clip_range")
        ncdf_varget, id, id1, crange
        id1 = ncdf_varid(id,"DZ")
        ncdf_varget, id, id1, dz
        dz = dz*0.01
        id1 = ncdf_varid(id,"VE")
        ncdf_varget, id, id1, ve
        ve = ve*0.01
        id1 = ncdf_varid(id,"SW")
        ncdf_varget, id, id1, sw
        sw = sw*0.01
        id1 = ncdf_varid(id,"time_offset")
        ncdf_varget, id, id1, offset

        NCDF_DIMINQ, id, 0, name, nazi

        x = fltarr(960,nazi,20)
        y = fltarr(960,nazi,20)
        z = fltarr(960,nazi,20)
        dbz = fltarr(960,nazi,20)
        vel = fltarr(960,nazi,20)
        spw = fltarr(960,nazi,20)
        range = fltarr(960,nazi,20)


        pi = 4.0*atan(1.0)

        a = 6378.137                            		;Earth's equatorial radius, km
        b = 6356.7523                           		;Earth's polar radius, km

        phi = (lat*pi)/180.                 			;radar latitude, radians

        top = (a^2*cos(phi))^2 + (b^2*sin(phi))^2
        bot = (a*cos(phi))^2 + (b*sin(phi))^2

        er = sqrt(top/bot)                      		;Earth radius at radar latitude, km (see Wikipedia "Earth Radius")

        eer = (4./3.)*er                        		;Effective earth radius (corrects for atmospheric refraction), km

	range(0,*,k) = (-1.*range2fc)/1000.			;km	

	for j=0,nazi-1 do begin
	last = ((crange(j) - (-1.*range2fc))/cellsp) + 1
	for i=1,960-1 do begin

	if (i le last) then begin
	
	range(i,j,k) = range(i-1,j,k) + cellsp/1000.		;Range to center of pulse volume, km	
        fi = (ele(j)*pi)/180.                           	;radians

        top = range(i,j,k)*cos(fi)
        bot = eer + range(i,j,k)*sin(fi)

        s = eer*atan(top/bot)                           	;km

;       An azimuth of 0 deg points to true north while 90 deg points east (i.e. clockwise rotation).    
        theta = (azi(j)*pi)/180.                        	;radians

        x(i,j,k) = s*sin(theta)               			;km, radar relative
        y(i,j,k) = s*cos(theta)                            	;km, radar relative
	dbz(i,j,k) = dz(i,j)					;reflectivity, dBZ
	vel(i,j,k) = ve(i,j)					;radial velocity, m/s
	spw(i,j,k) = sw(i,j)					;spectral width, m/s

	lam = atan(top/bot)					;radians
	
	psi = pi - (lam + pi/2.)				;radians

	hyp = eer/cos(lam)					;km

	p1 = (range(i,j,k)*sin(fi))/sin(psi)			;km
	p2 = hyp - eer						;km

	z(i,j,k) = p1 + p2 + (alt/1000.)			;km, Earth relative

	foot = sin((0.95*pi)/180.)*range(i,j,k)			;km, beam footprint


	endif else begin

	range(i,j,k) = -327.680
	x(i,j,k) = -327.680
	y(i,j,k) = -327.680
	z(i,j,k) = -327.680
	dbz(i,j,k) = -327.680
	vel(i,j,k) = -327.680
	spw(i,j,k) = -327.680
	
	endelse

	endfor	
	endfor	

;	writeu,3,x(*,*,k),y(*,*,k),z(*,*,k),dbz(*,*,k)

	ze = 10^(dbz/10.)               ;Z units (for interpolation because its linear)

	if (k eq 0 or k eq 4 or k eq 8 or k eq 11 or k eq 14 or $ 
k eq 16 or k eq 17 or k eq 18 or k eq 19) then begin

	cnt = cnt + 1 

	hold = GRIDXY(x(*,*,k),y(*,*,k),ze(*,*,k))

	car_r(*,*,cnt) = hold

	hold = GRIDXY(x(*,*,k),y(*,*,k),z(*,*,k))
	
	car_z(*,*,cnt) = hold

;	hold = GRIDXY(x(*,*,k),y(*,*,k),vel(*,*,k))
;	car_v(*,*,k) = hold

	endif

	close,2

	print,'FINISHED SCAN:',k+1
	ENDFOR

	openw,8,'ferret_input'
	writeu,8,car_z,car_r

;	gx = fltarr(961,20)
	
;        xin = -240.
;	for k = 0,19 do begin
;        for xx = 0,960 do begin         ;making grid points in x
	
;        gx[xx,k] =  xin + xx*0.5

;        endfor
;	endfor

;	Put in perturbation to get triangulate to work
;	seed = 10001L
;	random = randomu(seed,961,20)
;	gx = gx + random

;	for j = 0,960 do begin
;	hold = GRIDZ(gx,car_z[*,j,*],car_r[*,j,*])

;	print,'FINISHED PLANE:',j+1
;	data_3d[*,j,*] = hold		;3-D gridded radar field
;	endfor
	
;	openw,5,'data_3d.dat'
;	writeu,5,data_3d

END