d2r = pi/180;
%Storm lat/lon at a central time...
slat = 21.680;slon = 360 - 58.303;

%lower_fuselage;	%subroutine to display NOAA LF data.

%Track-relative grid spacing
delx = 0.5;     %km
dely = 0.5;     %km
delz = 0.030;   %km

rx = -4:delx:4;
%ry = -1:dely:10;
rz = delz:delz:2.5;

nx = numel(rx);
%ny = numel(ry);
nz = numel(rz);

naz1 = zeros(nx,nz);
naz2 = zeros(nx,nz);
naz3 = zeros(nx,nz);
%nel1 = zeros(nx,ny,nz);
%nel2 = zeros(nx,ny,nz);
%nel3 = zeros(nx,ny,nz);

y = 20;

ix = 0;
for x = rx
    ix = ix + 1;
    iy = 0;
    x
    tic
%    for y = ry
%        iy = iy + 1;
        iz = 0;
        for z = rz
            iz = iz + 1;
 
            
fid = fopen('test.txt');
g=textscan(fid,'%s');

for a=1:4
path = fullfile(g{1}{a});
ncid = netcdf.open(path,'NC_NOWRITE');
[ndims,nvars,natts,unlimdimID]= netcdf.inq(ncid);

%Variable indices...
%0 = time, 1 = incidence angle, 2 = a/c lon, 3 = a/c lat, 4 = lon @ range gate 100, 5 = lat @ range gate 100,
%6 = a/c altitude (meters), 7 = azimuth w.r.t. true N, 8 = track (deg), 9 = yaw (deg), 10 = pitch (deg), 11 = roll (deg),
%12 = heading (deg)

%[dimname, dimlen] = netcdf.inqDim(ncid,0);
%[varname, vartype, dimids, natts] = netcdf.inqVar(ncid,0);

unixt = netcdf.getVar(ncid,0);		%Unix time stamp (seconds since Jan 1, 1970 00:00 UTC
test = datenum([1970 1 1 0 0 unixt(1,1)]);
time = datestr(test);

iang = netcdf.getVar(ncid,1);
alon = netcdf.getVar(ncid,2);
alat = netcdf.getVar(ncid,3);
rlon = netcdf.getVar(ncid,4);
rlat = netcdf.getVar(ncid,5);
ralt = netcdf.getVar(ncid,6);
azi = netcdf.getVar(ncid,7);
track = netcdf.getVar(ncid,8);
pitch = netcdf.getVar(ncid,10);
roll = netcdf.getVar(ncid,11);
head = netcdf.getVar(ncid,12);

drift = track - head;

dis = zeros(1,size(track,2));

%incidence;		%subroutine to compute nominal incidence angle
%end
%tau = mean(vals);	%mean nominal incidence angle

%nominal incidence angles below (variable includes 1.5 degree pitch offset for antenna mount)

%KU band...4 beams
%*-0-.nc files, outer beam
%tau = 48.60	%degrees
%*-1-.nc files
%tau = 40.50	%degrees
%*-2-.nc files
%tau = 36.20	%degrees
%*-3-.nc files, inner beam
%tau = 32.30	%degrees

%C band...4 beams
if (a == 1)
    %*-0-.nc files, outer beam
    tau = 50.00;     %degrees
elseif (a == 2)
    %*-1-.nc files
    tau = 44.10;     %degrees
elseif (a == 3)
    %*-2-.nc files
    tau = 39.10;     %degrees
elseif (a == 4)
    %*-3-.nc files, inner beam
    tau = 34.50;     %degrees
end


%Need to find the range gate that intersects the surface.  Without pitch and roll, this can be calculated
%based on the nominal incidence angle and the altitude of the aircraft.  See how bad this assumption will be.
%Just use instantaneous inci angle and see how much variability in sfc there is, then use min value.

rsfc = ralt./cos(iang*d2r);
sfc_rg = round(mean(mean(rsfc))/30);            %mean surface range gate, 30 m gate spacing	

range = ((1:1:sfc_rg)*30)/1000;                 %km
range = repmat(range,size(head,1)*size(head,2),1);
range = reshape(range,size(head,1),size(head,2),sfc_rg);

azimuth = repmat(azi,[1 1 sfc_rg]);

ralt = ralt/1000;       %km

%compute x,y sfc position (relative to aircraft) of each binned measurement, meters
%radius = ralt.*tan(iang.*d2r);	%radius of surface footprint, m
%x = radius.*cos(azi);
%y = radius.*sin(azi);
%flat = alat + y/(111*1000);				%degrees
%flon = (alon+360.) + x./(111*1000*cos(flat.*d2r));	%degrees

%reply = input('Do you want a storm-centered grid? Y/N: ', 's');

reply = 'N';
stormcentered = strcmp(reply,'Y');

%mapearth;

if (stormcentered == 1)

disp('Storm Centered')
mapearth;	%map data from aicraft to earth coordinates
%Below is for storm-centered grid, slat/slon are center coords 
%-------------------------------------------------------------
%mxlat = max(max(ilat(:,:,108)));mxlon = max(max(ilon(:,:,108)));
%slat = 21.680;%py = (mxlat - slat)*111;
%slon = 360 - 58.303;%px = (mxlon - slon)*111*cos(slat*d2r);
%Making grid
resxy = 0.20;                           %km
%xc = [-round(max(py,px)):resxy:round(max(py,px))];yc = xc;     %km
xc = -2:resxy:15;yc = -40:resxy:10;%km
resz = 0.030;                           %km
zc = resz:resz:2.40;
[xc,yc,zc]=meshgrid(xc,yc,zc);
%Cartesian storm-centered radar coordinates
xr = (ilon - slon)*111*cos(slat*d2r);   %km
yr = (ilat - slat)*111;                 %km
zr = ihgt;                              %km
%-------------------------------------------------------------

else

%Below is for aircraft-pass-centered grid
%disp('Aircraft Track Centered')

maptrack;	%map data from aicraft to track coordinates

dis(1) = 0.;		%start of pass
for i=1:size(range,2)-1
con = 111*cos(alat(1,i)*d2r);
dy = (alat(1,i+1) - alat(1,i))*111;			%km
dx = (alon(1,i+1)*con - alon(1,i)*con);			%km
dis(i+1) = (dx^2 + dy^2)^0.5 + dis(i);			%km
end

%dis = repmat(dis,[32 1]);

%Making grid
resxy = 0.20;	%km
xc = -5:resxy:5;yc = -5:resxy:50;	%km
resz = 0.030;                           %km
zc = resz:resz:2.40;
[xc,yc,zc] = meshgrid(xc,yc,zc);
%Cartesian track-centered radar coordinates
xr = Xt;                                %km
yr = bsxfun(@plus,Yt,dis);              %km
zr = ihgt;                              %km
%-------------------------------------------------------------

end

         
%Keeping data only in grid cell of interest
mask = xr >= (x-delx/2) & xr < (x+delx/2) & yr >= (y-dely/2) & yr < (y+dely/2) & zr >= (z-delz/2) & zr < (z+delz/2);

xr = xr(mask);
yr = yr(mask);
zr = zr(mask);
Yt = Yt(mask);

%xr(~mask) = NaN;
%yr(~mask) = NaN;
%zr(~mask) = NaN;

%Xe(~mask) = NaN;
%Ye(~mask) = NaN;
%Ze(~mask) = NaN;

%ind = find(mask == 1);
%x3 = xr(ind);
%y3 = yr(ind);
%z3 = zr(ind);
%[i,j,k]=ind2sub([32,399,115],ind);

%xr = reshape(xr,size(xr,1)*size(xr,2)*size(xr,3),1);
%yr = reshape(yr,size(xr,1)*size(xr,2)*size(xr,3),1);
%zr = reshape(zr,size(xr,1)*size(xr,2)*size(xr,3),1);
%Yt = reshape(Yt,size(xr,1)*size(xr,2)*size(xr,3),1);

%azimuth3 = reshape(azimuth(i,:,:),size(xr,2)*size(xr,3),1);
    
%xe = reshape(Xe(i,:,:),size(xr,2)*size(xr,3),1);
%ye = reshape(Ye(i,:,:),size(xr,2)*size(xr,3),1);
    
%Earth relative spherical coords
%ran = sqrt(xe.^2 + ye.^2 + z3.^2);  %range (km)
%azi = atan(xe./ye)*(1/d2r);         %azimuth (degrees)
%ele = asin(z3./ran)*(1/d2r);        %elevation (degrees)

%Track relative spherical coords
    
ran = sqrt(xr.^2 + Yt.^2 + zr.^2);    %range (km)
ele = acos(zr./ran)*(1/d2r);          %elevation (degrees)
    
yp = Yt > 0.0;
foreb = 90 + atan(Yt(yp)./xr(yp))*(1/d2r);   %fore azimuth (degrees)   
    
yn = Yt < 0.0;
aftb = atan(xr(yn)./Yt(yn))*(1/d2r);         %aft azimuth(degrees)
    

if (a == 1)
    cxr = xr;
    cyr = yr;
    czr = zr;
    cazi = cat(1,foreb,aftb);
    cele = ele;
else
    cxr = cat(1,cxr,xr);
    cyr = cat(1,cyr,yr);
    czr = cat(1,czr,zr);
    cazi = cat(1,cazi,cat(1,foreb,aftb));
    cele = cat(1,cele,ele);
end


netcdf.close(ncid);
end

%scatter3(cxr,cyr,czr,50,cazi,'filled');

save = 1000;
%Finding optimal azimuth diversity
for b = 1:1:size(cazi,1)
    tmp = cazi;     %save original azimuths
    num = tmp(b);   %select one azimuth to examine
    tmp(b) = [];    %remove that azimuth from original array so we don't subtract it from itself
    [bazi,ind] = min(abs(tmp - num - 90));  %subtract selected azimuth from rest of azimuths in array relative to optimal (90), then minimize it.
    val = tmp(ind);
    if (bazi < save)
        save = bazi;
        saz1 = num;
        saz2 = val;
    end
end

if (numel(cazi) == 0)
    naz1(ix,iz) = NaN;
    naz2(ix,iz) = NaN;
    naz3(ix,iz) = NaN;
%    nel1(ix,iz) = NaN;
%    nel2(ix,iz) = NaN;
%    nel3(ix,iz) = NaN;
else
    naz1(ix,iz) = save;
    naz2(ix,iz) = saz1;
    naz3(ix,iz) = saz2;
    
%    naz1(ix,iy,iz) = nanmean(cazi);                         %mean azimuth look in cell (deg)
%    naz2(ix,iy,iz) = nanstd(cazi);                          %std dev of azimuth look in cell (deg)
%    naz3(ix,iy,iz) = abs(max(cazi) - min(cazi));            %range of azimuth look in cell (deg)
%    nel1(ix,iy,iz) = nanmean(cele);                         %mean elevation look in cell (deg)
%    nel2(ix,iy,iz) = nanstd(cele);                          %std dev of elevation look in cell (deg)
%    nel3(ix,iy,iz) = abs(max(cele) - min(cele));            %range of elevation look in cell (deg)
end

%scatter3(nxr,nyr,nzr,50,naz,'filled');
%hold on

%count(ix,iy,iz) = sum(sum(sum(mask)));
%        end
%    end
%end

fclose(fid);

        end
%    end
    toc
end

%contourf(rx,ry,naz3(:,:,1)')
%xlabel('Across Track');ylabel('Along Track');
%return;



return;

xlim([(x-0.1) (x+0.1)]);
ylim([(y-0.1) (y+0.1)]);
zlim([(z-0.015) (z+0.015)]);



return;

%ind = find(mask == 1);
%x3 = xr(ind);
%y3 = yr(ind);
%z3 = zr(ind);
%[i,j,k]=ind2sub([32,399,115],ind);
%count(ix,iy,iz) = sum(sum(sum(mask)));









%Get reflectivity and Doppler data

%path = fullfile('/mercury/guimond/GSFC/IWRAP/20030912/volume_reflectivity/',z{1}{16})
path = fullfile('/mercury/guimond/GSFC/IWRAP/20030912/volume_reflectivity/',z{1}{a})
ncid = netcdf.open(path,'NC_NOWRITE');

%Variable indices...
%0 = reflectivity, 1 = reflectivity noise removed
refl = netcdf.getVar(ncid,1); 	%units m^6/m^3
missing = refl == 0.0;         	%if inequality is true, missing is 1; false, missing is 0.
refl(missing) = NaN;     	%assign NaN to places in original array where missing is 1, otherwise use original data

dbz = 10.*log10(abs(refl)) + 180;	%units of dBZ

dbz = shiftdim(dbz(3:sfc_rg,:,:),1);	%NOTE: DO NOT use reshape here, will mess up array.

ndbz = double(reshape(dbz,size(xr,1)*size(xr,2)*size(xr,3),1));
nxr = double(reshape(xr,size(xr,1)*size(xr,2)*size(xr,3),1));
nyr = double(reshape(yr,size(yr,1)*size(yr,2)*size(yr,3),1));
nzr = double(reshape(zr,size(zr,1)*size(zr,2)*size(zr,3),1));

F = TriScatteredInterp(nxr,nyr,nzr,ndbz);
dbzg = F(xc,yc,zc);
dbzg = smooth3(dbzg);
%test = squeeze(dbzg(:,50,:));

%dbzg = griddata3(double(xr),double(yr),double(zr),double(dbz),double(xc),double(yc),double(zc),'linear');

%set(gcf,'PaperPositionMode','auto');	sets printed figure to be same size as window figure.
%set(gcf, 'PaperUnits', 'inches');  
%set(gcf, 'PaperSize', [8.5 11]);   
%set(gca,'XTick',[-5 -4 -3 -2 -1 0 1 2 3 4 5])

return;

d = 12;
%for d=1:32
clearvars dis height
dis(1) = 0.;		%start of pass
for i=1:398
con = 111*cos(rlat(d,i)*d2r);
dy = (rlat(d,i+1) - rlat(d,i))*111;			%km
dx = (rlon(d,i+1)*con - rlon(d,i)*con);			%km
dis(i+1) = (dx^2 + dy^2)^0.5 + dis(i);			%km
end



s = reshape(dbz(:,d,:),200,399);

%if (a == 1)
%   pass = s;
%elseif (a == 5 || a == 9 || a == 13 || a == 17)
%   pass = cat(2,pass,s);
%end;

pass = s;

%end

sfc_rg = 100;

%first useful range gate = 3, surface range gate = 81
%height = (([3:1:81]-3)*30)/1000;	%km
height = (([3:1:sfc_rg]-3)*30)/1000;	%km
[dis,height] = meshgrid(dis,height);

pass = flipdim(pass(3:sfc_rg,:),1);

v = 30:2:52;

%contourf(flipud(s))
%set(gca,'YDir','rev') 

%%contourf(dis,height,pass,v)          

%set(gca,'YDir','rev') 
%%colorbar;caxis([30 52]);
%%ylim([0 ((81-3)*30/1000)])
%get(gcf,'Position')
%%set(gcf,'Position',[0 0 1200 400])
%%pbaspect([4 1 1])
%daspect([1 1 1])
%%xlabel('Along Track (km)');
%%ylabel('Height (km)');
%%title('IWRAP Isabel Eyewall 9/12/2003 Raw Data');

%print -depsc raw

avgf = fspecial('average');
new = imfilter(pass,avgf);

%Fixing boundaries
%new(1,:) = pass(1,:);
%new(79,:) = pass(79,:);
%new(:,1) = pass(:,1);
%new(:,399) = pass(:,399);

contourf(dis,height,new,v)          
colorbar;caxis([30 52]);
ylim([0 ((sfc_rg-3)*30/1000)])
set(gcf,'Position',[0 0 1200 400])
pbaspect([4 1 1])

xlabel('Along Track (km)');
ylabel('Height (km)');
title(['IWRAP (Azi pos. ' int2str(d) ' ) Isabel Eyewall 9/12/2003 3 point Running Mean Filter']);


%fname{d} = ['filter' int2str(d)];


%print('-dpng',fname{d})
%end

%missing = new == -9999;         %if inequality is true, missing is 1; false, missing is 0.
%new(missing) = NaN;             %assign NaN to places in original array where missing is 1, otherwise use original data
%lh = new(:,:,:,:,2);
