%This program uses the newest IWRAP lev 1b files (2/4/2011)

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.

fid = fopen('passes_cband_lev1b.txt');
g=textscan(fid,'%s');

%fid = fopen('/mercury/guimond/GSFC/IWRAP/20030912/volume_reflectivity/passes_cband.txt');
%z=textscan(fid,'%s');

cnt = 0;

gate = 30;      %gate spacing, m

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


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

% Get name of global attribute (nominal incidence angle in degrees)
gattname = netcdf.inqAttName(ncid,netcdf.getConstant('NC_GLOBAL'),9);

% Get value of global attribute (nominal incidence angle in degrees)
tau = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),gattname)

%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)]);
%time = datestr(test);

unixt = netcdf.getVar(ncid,0); %Unix time stamp (seconds since Jan 1, 1970 00:00 UTC
    yy = unixt*0+1970;
    mm = unixt*0+1;
    dd = unixt*0+1;
    mn = unixt*0;
    ss = unixt*0;
    time = datenum([yy mm dd mn ss unixt]);
%    dtg = datestr(time);    

azi = netcdf.getVar(ncid,1);
scan = netcdf.getVar(ncid,3);
nav = netcdf.getVar(ncid,4);
alongtrack = netcdf.getVar(ncid,5);
Ze = netcdf.getVar(ncid,6);
rho = netcdf.getVar(ncid,7);
vel = netcdf.getVar(ncid,8);
%nrcs = netcdf.getVar(ncid,9);
%mZe = netcdf.getVar(ncid,10);
svel = netcdf.getVar(ncid,11);
thetap = netcdf.getVar(ncid,12);
%thetam = netcdf.getVar(ncid,13);
%rdist = netcdf.getVar(ncid,14);
%flags = netcdf.getVar(ncid,15);
mask = netcdf.getVar(ncid,16);
vcorr = netcdf.getVar(ncid,17);
%trcs = netcdf.getVar(ncid,18);


nz = size(Ze,1);    %number of range gates
pad = 200 - nz;

%Below is for Sep 12, 2003 1854 UTC file...
%time = time(2:length(scan));
%azi = azi(2:length(scan));
%Ze = Ze(2:length(scan),:);
%vel = vel(:,2:length(scan));
%svel = svel(2:length(scan));
%mask = mask(2:length(scan));
%vcorr = vcorr(2:length(scan));
%scan = scan(2:length(scan));


ii = 1;     %first point is start of azimuthal scan
for i=1:31
    if (scan(i)~= scan(i+1)) 
        ii = i+1;
    end
end


ff = length(scan);  %lst point is end of azimuthal scan
for i = ff:-1:ff-30
    if (scan(i) ~= scan(i-1))
        ff = i-1;
    end
end

pp = (ff-ii)+1;

%dtg = dtg(ii:ff,:);dtg = reshape(dtg,32,pp/32,20);

alongtrack = alongtrack(ii:ff);
alongtrack = reshape(alongtrack,32,pp/32);

time = time(ii:ff);
time = double(reshape(time,size(time,1)*size(time,2),1));

rho = shiftdim(rho(:,ii:ff),1);
rho = reshape(rho,32,pp/32,nz);

%For Isabel only ?
iwrap_reflectivity;         %get calibrated reflectivity from older files

%removing duplicates, keeping first instance
[zzz,ind] = unique(ctime,'first');
ctime = ctime(sort(ind));   %reverse unique's sorting
cdbz = cdbz(sort(ind),:);

clear dbzn

%Interpolating calibrated reflectivity times to new file times
for i=1:nz
    dbzn(:,i) = interp1(ctime,cdbz(:,i),time,'linear');
end

mask = shiftdim(mask(:,ii:ff),1);
mask = reshape(mask,32,pp/32,nz);

dbzn = reshape(dbzn,32,pp/32,nz);  %calibrated reflectivity, Z
dbzn = padarray(dbzn,[0 0 pad],'post');


%Cleaning up reflectivity...can also just cut 0.5 km above ocean surface
%ind = find(mask ~= 1);
%dbzn(ind) = NaN;
%ind = find(10*log10(dbzn) > 55.0);
%dbzn(ind) = NaN;


%Reflectivity processing...
%missing = Ze <= 0.0;      %if inequality is true, missing is 1; false,
%missing is 0.
%Ze = shiftdim(Ze(:,ii:ff),1);
%Ze = reshape(Ze,32,pp/32,nz);
%missing = Ze == 0.0;
%Ze(missing) = NaN;     
%Ze = abs(Ze);
%dbz = 10*log10(Ze);
%dbz = dbz(ii:ff,:);
%dbz = shiftdim(dbz(:,ii:ff),1);
%dbz = reshape(dbz,32,pp/32,nz);

%avgf = fspecial('average');
%pass = squeeze(dbzn(15,:,:));
%pass = 10*log10(pass);          %dBZ
%sdbz = imfilter(pass,avgf);
%v = 30:1:55;
%contourf(sdbz',v)
%set(gca,'YDir','reverse')

vel = shiftdim(vel(:,ii:ff),1);
vel = reshape(vel,32,pp/32,nz);

vcorr = vcorr(ii:ff);
vcorr = reshape(vcorr,32,pp/32);

cvel = bsxfun(@minus,vel,vcorr);    %Doppler velocity with aircraft motion removed
cvel = padarray(cvel,[0 0 pad],'post');

%cvel(isnan(dbzn)) = NaN;

%Computing aircraft motion correction on my own...
%Taut = asin(bsxfun(@rdivide,Yt,reshape(range,1,1,size(range,2))));
%Note:  ground speed not in nav file.
%mycorr = bsxfun(@times,sin(Taut),nav(:,:,13)) +
%bsxfun(@times,Zt,nav(:,:,7));

%svel = svel(ii:ff);
%svel = reshape(svel,32,pp/32);
%fix2 = svel - vcorr;

nav = shiftdim(nav,1);
nav = nav(ii:ff,:);
nav = reshape(nav,32,pp/32,34);

azi = azi(ii:ff);
azi = reshape(azi,32,pp/32);

range = ([1:1:nz]*gate)/1000;     %km
range = repmat(range,pp,1);
range = reshape(range,32,pp/32,nz);

ralt = nav(:,:,15)/1000;        %km (radar altitude)
roll = nav(:,:,11);             %deg
pitch = nav(:,:,10);            %deg
drift = nav(:,:,8);             %deg

alat = nav(:,:,3);              %deg
alon = nav(:,:,4)+360;          %deg

alat = mean(alat,1);
alon = mean(alon,1);

head = nav(:,:,9);              %deg

%%figure(5)
%%hold on
%plot(alon,alat,'k','linewidth',4)

%%ilon = alon <= max(max(lon)) & alon >= min(min(lon));

%%alon = alon(ilon);
%%alat = alat(ilon);

%%xa = [alon(1) alon(size(alon,2))];
%%ya = [alat(1) alat(size(alat,2))];

%%[xaf,yaf] = ds2nfu(xa,ya);
%%annotation('arrow',xaf,yaf,'LineWidth',3)

%print('-dpdf',[datestr(time(1)) 'LFscan'])
%close(5)

%%maptrack;
%%Xt = padarray(Xt,[0 0 pad],'post');
%%ihgt = padarray(ihgt,[0 0 pad],'post');

%assuming alongtrack distance not a function of azimuth, error of ~ 5 m.
tdis = alongtrack(1,:) - alongtrack(1,1);   %along-track distance, km 

ax = bsxfun(@times,sin(head*d2r),tdis);
ay = bsxfun(@times,cos(head*d2r),tdis);

%dis(1) = 0.;		%start of pass
%for i=1:599-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

%Track-relative coordinates centered on radar
%%yri = bsxfun(@plus,Yt,tdis);	%km
%%yri = padarray(yri,[0 0 pad],'post');

mapearth;
Xe = padarray(Xe,[0 0 pad],'post');     %km
Ye = padarray(Ye,[0 0 pad],'post');     %km
Ze = padarray(Ze,[0 0 pad],'post');     %km
ihgt = padarray(ihgt,[0 0 pad],'post'); %km

aXe = bsxfun(@plus,Xe,ax);       %km
aYe = bsxfun(@plus,Ye,ay);       %km

range = padarray(range,[0 0 pad],'post');

if (b == 1)
    Vr = cvel;
    dbz = dbzn;
    xr = Xe;
    yr = Ye;
    zr = Ze;
    xxe = aXe;
    yye = aYe;
    zze = ihgt;
    nrange = range;
elseif (b == 2)
    Vr = cat(4,Vr,cvel);
    dbz = cat(4,dbz,dbzn);
    xr = cat(4,xr,Xe);
    yr = cat(4,yr,Ye);
    zr = cat(4,zr,Ze);
    xxe = cat(4,xxe,aXe);
    yye = cat(4,yye,aYe);
    zze = cat(4,zze,ihgt);
    nrange = cat(4,nrange,range);
end

end


delx = 0.50;        %km
dely = 0.50;        %km
delz = gate/1000;   %km

%Track-relative grid, km
rx = -4:delx:4;nx = numel(rx);                  %across-track
ry = 0:dely:70;ny = numel(ry);                  %along-track
rz = 0.5:delz:2.5;nz = numel(rz);
%rz = 1:1:10;nz = numel(rz);

nhead = interp1(tdis,head(1,:),ry,'linear');

rz = shiftdim(rz,-1);
cz = repmat(rz,[ny nx 1]);

rx = repmat(rx,[ny 1 nz]);

rhx = ry.*sin(nhead*d2r);rhx = shiftdim(rhx);
nrhx = repmat(rhx,[1 nx nz]);
rhy = ry.*cos(nhead*d2r);rhy = shiftdim(rhy);
nrhy = repmat(rhy,[1 nx nz]);

nhead = shiftdim(nhead);
nhead = repmat(nhead,[1 nx nz]);

cx = nrhx + rx.*cos(nhead*d2r);
cy = nrhy + rx.*sin(nhead*d2r)*-1;

nrange = nrange*1000;   %range, m

%Looping over grid cells...using track-relative coords projected onto earth
%coords

ur = zeros(ny,nx,nz)+NaN;
vr = zeros(ny,nx,nz)+NaN;
Wr = zeros(ny,nx,nz)+NaN;
gdbz = zeros(ny,nx,nz)+NaN;

rsq = zeros(ny,nx,nz)+NaN;
count = zeros(ny,nx,nz)+NaN;

dr = 1.0;     %search radius, km

%zz = [1 5 10];

%zz = 1:1:nz;

%zz = 1;

zz = squeeze(squeeze(rz));

for iy = 51:ny
    iy  %#ok<NOPTS>
    tic
    for ix = 1:nx
%        iz = 0;
%        iz = 1;
        for iz = 1:nz
%        iz = iz + 1;

%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);

x = cx(iy,ix,iz);
y = cy(iy,ix,iz);
%z = cz(iy,ix,iz);
z = zz(iz);

radius = sqrt((xxe - x).^2 + (yye - y).^2);

mask = radius <= dr & zze >= (z-delz/2) & zr < (z+delz/2);

%%xr2 = xr(mask);
%%yr2 = yr(mask);
%%zr2 = zr(mask);
%%yy2 = yy(mask);

Xe = xr(mask)*1000;     %m
Ye = yr(mask)*1000;     %m
Ze = zr(mask)*1000;     %m

Y = Vr(mask).*nrange(mask);	%radial velocity observations times range, m^2/s
X = [Xe,Ye,Ze];             %Cartesian Earth coordinates of radar observations, m

B = X\Y;                    %Solving least squares problem, B = [u, v, W], m/s

if (numel(Y) < 3)
ur(iy,ix,iz) = NaN;
vr(iy,ix,iz) = NaN;
Wr(iy,ix,iz) = NaN;
gdbz(iy,ix,iz) = NaN;

rsq(iy,ix,iz) = NaN;
count(iy,ix,iz) = NaN;
else
%weights = abs(radius(mask) - dr)/dr;       %linear distance weights
%weights = exp(-1*(radius(mask).^2/(2*dr^2)));     %Barnes (Gaussian) distance weight   
%B = lscov(X,Y,weights);     %Solving weighted least squares problem
    
ur(iy,ix,iz) = B(1);
vr(iy,ix,iz) = B(2);
Wr(iy,ix,iz) = B(3);
gdbz(iy,ix,iz) = 10*log10(nanmean(dbz(mask)));


sstot = nansum((Y - nanmean(Y)).^2);
sserr = nansum((Y - X*B).^2);

rsq(iy,ix,iz) = 1 - (sserr./sstot);     %explained variance
count(iy,ix,iz) = numel(Y);             %number of obs in grid cell
end

%numel(Y)

%rank(Y)
%cond(Y)
        end
    end
    toc
end


return;





%str = fprintf('location fore (degrees,xtrack(km)): %f, %f\n',azi(locf),mxp);
%str = fprintf('location aft (degrees,xtrack(km)): %f, %f\n',azi(loca),mxn);

forex(c) = mxp;
aftx(c) = mxn;

if (c == 1)
    xf = Xt(locf,c,:);xa = Xt(loca,c,:);
    yf = Yt(locf,c,:);ya = Yt(loca,c,:);
    yyf = yr(locf,c,:);yya = yr(loca,c,:);
    zf = Zt(locf,c,:);za = Zt(loca,c,:);
    zzf = ihgt(locf,c,:);zza = ihgt(loca,c,:);    
    reflf = dbzn(locf,c,:);refla = dbzn(loca,c,:);
    dvelf = cvel(locf,c,:);dvela = cvel(loca,c,:);
    rhof = rho(locf,c,:);rhoa = rho(loca,c,:);
%    maskf = mask(locf,c,:);maska = mask(loca,c,:);
else
    xf = cat(3,xf,Xt(locf,c,:));xa = cat(3,xa,Xt(loca,c,:));
    yf = cat(3,yf,Yt(locf,c,:));ya = cat(3,ya,Yt(loca,c,:));
    yyf = cat(3,yyf,yr(locf,c,:));yya = cat(3,yya,yr(loca,c,:));
    zf = cat(3,zf,Zt(locf,c,:));za = cat(3,za,Zt(loca,c,:));
    zzf = cat(3,zzf,ihgt(locf,c,:));zza = cat(3,zza,ihgt(loca,c,:));
    reflf = cat(3,reflf,dbzn(locf,c,:));refla = cat(3,refla,dbzn(loca,c,:));
    dvelf = cat(3,dvelf,cvel(locf,c,:));dvela = cat(3,dvela,cvel(loca,c,:));
    rhof = cat(3,rhof,rho(locf,c,:));rhoa = cat(3,rhoa,rho(loca,c,:));
%    maskf = cat(3,maskf,mask(locf,c,:));maska = cat(3,maska,mask(loca,c,:));
end




%Making grid (km in x/y/z)
resxy = 0.20;          
endt = ceil(max(tdis)) + 2;
xc1 = 0;yc1 = [0:resxy:endt];yc1 = double(yc1);	
resz = 0.030;               
zc1 = [resz:resz:2.40];     
[xc,yc,zc]=meshgrid(xc1,yc1,zc1);

xf = double(reshape(xf,size(xf,1)*size(xf,2)*size(xf,3),1));
xa = double(reshape(xa,size(xf,1)*size(xf,2)*size(xf,3),1));
yf = double(reshape(yf,size(xf,1)*size(xf,2)*size(xf,3),1));
ya = double(reshape(ya,size(xf,1)*size(xf,2)*size(xf,3),1));
yyf = double(reshape(yyf,size(xf,1)*size(xf,2)*size(xf,3),1));
yya = double(reshape(yya,size(xf,1)*size(xf,2)*size(xf,3),1));
zf = double(reshape(zf,size(xf,1)*size(xf,2)*size(xf,3),1));
za = double(reshape(za,size(xf,1)*size(xf,2)*size(xf,3),1));
zzf = double(reshape(zzf,size(xf,1)*size(xf,2)*size(xf,3),1));
zza = double(reshape(zza,size(xf,1)*size(xf,2)*size(xf,3),1));
reflf = double(reshape(reflf,size(xf,1)*size(xf,2)*size(xf,3),1));
refla = double(reshape(refla,size(xf,1)*size(xf,2)*size(xf,3),1));
dvelf = double(reshape(dvelf,size(xf,1)*size(xf,2)*size(xf,3),1));
dvela = double(reshape(dvela,size(xf,1)*size(xf,2)*size(xf,3),1));
rhof = double(reshape(rhof,size(xf,1)*size(xf,2)*size(xf,3),1));
rhoa = double(reshape(rhoa,size(xf,1)*size(xf,2)*size(xf,3),1));
%maskf = double(reshape(maskf,size(xf,1)*size(xf,2)*size(xf,3),1));
%maska = double(reshape(maska,size(xf,1)*size(xf,2)*size(xf,3),1));

rf = sqrt(yf.^2 + zf.^2);       %km
ra = sqrt(ya.^2 + za.^2);       %km

%Gridding fore beam data...
F = TriScatteredInterp(xf,yyf,zzf,reflf);
dbzfg = F(xc,yc,zc);
dbzfg = 10*log10(dbzfg);        %dBZ
dbzfg = squeeze(dbzfg);

F = TriScatteredInterp(xf,yyf,zzf,dvelf);
dopfg = F(xc,yc,zc);            %m/s

F = TriScatteredInterp(xf,yyf,zzf,zf);
zfg = F(xc,yc,zc)*1000;         %m

F = TriScatteredInterp(xf,yyf,zzf,yf);
yfg = F(xc,yc,zc)*1000;         %m

F = TriScatteredInterp(xf,yyf,zzf,rf);
rfg = F(xc,yc,zc)*1000;         %m

F = TriScatteredInterp(xf,yyf,zzf,rhof);
rhofg = F(xc,yc,zc);            %dimensionless
rhofg = squeeze(rhofg);

%F = TriScatteredInterp(xf,yyf,zzf,maskf,'nearest');
%mfg = F(xc,yc,zc);
%mfg = squeeze(mfg);

%Gridding aft beam data...
F = TriScatteredInterp(xa,yya,zza,refla);
dbzag = F(xc,yc,zc);
dbzag = 10*log10(dbzag);        %dBZ
dbzag = squeeze(dbzag);

F = TriScatteredInterp(xa,yya,zza,dvela);
dopag = F(xc,yc,zc);            %m/s

F = TriScatteredInterp(xa,yya,zza,za);
zag = F(xc,yc,zc)*1000;         %m

F = TriScatteredInterp(xa,yya,zza,ya);
yag = F(xc,yc,zc)*1000;         %m

F = TriScatteredInterp(xa,yya,zza,ra);
rag = F(xc,yc,zc)*1000;         %m

F = TriScatteredInterp(xa,yya,zza,rhoa);
rhoag = F(xc,yc,zc);            %dimensionless
rhoag = squeeze(rhoag);

%F = TriScatteredInterp(xa,yya,zza,maska,'nearest');
%mag = F(xc,yc,zc);
%mag = squeeze(mag);

ppcc = (rhofg + rhoag)/2;

vertical = (dopag.*rag.*yfg - dopfg.*rfg.*yag)./(yfg.*zag - yag.*zfg);     %m/s
vertical = squeeze(vertical);
%vertical = vertical + vt;                     %m/s

along = (dopag.*rag.*zfg - dopfg.*rfg.*zag)./(yag.*zfg - yfg.*zag);         %m/s
along = squeeze(along);

%Cleaning up data where there is contamination...
ind = find(zc1 < 0.50);
ind = max(ind)+1;
dbzfg = dbzfg(:,ind:size(zc1,2));
dbzag = dbzag(:,ind:size(zc1,2));
vertical = vertical(:,ind:size(zc1,2));
along = along(:,ind:size(zc1,2));
ppcc = ppcc(:,ind:size(zc1,2));
zc1 = zc1(:,ind:size(zc1,2));

dbz = (dbzfg + dbzag)/2;
fallspeeds;         %computing hydrometeor fallspeeds

w = vertical + vr;

%Masking data according to pulse-pair correlation coefficient
%(signal-to-noise)

ind = find(ppcc < 0.25 | dbz < 32);
dbz(ind) = NaN;
w(ind) = NaN;
along(ind) = NaN;

avgf = fspecial('average');
dbz2 = imfilter(dbz,avgf,'replicate');
ind = find(isnan(dbz2));
dbz2(ind) = dbz(ind);   %this step fixes boundary problem when filtering
w2 = imfilter(w,avgf,'replicate');
w2(ind) = w(ind);
along2 = imfilter(along,avgf,'replicate');
along2(ind) = along(ind);


%Plotting below...

endd = size(time,1);

%%figure(1)
%%set(gcf,'Position',[0 0 1200 300])
%%v = 30:1:52;
%%contourf(yc1,zc1,dbz2',v)
%%ylabel('Height (km)');
%%xlabel('Along Track (km)');
%%title(['Start:  ',datestr(time(1)),'  End:  ',datestr(time(endd))])
%%colorbar;

%%set(gcf,'PaperPositionMode','auto');	%sets printed figure to be same size as window figure.
%%set(gcf, 'PaperUnits', 'inches');  
%%set(gcf, 'PaperSize', [11 8.5]);   


%print('-dpdf',[datestr(time(1)) 'reflectivity'])
%close(1)
%print -dpng iwrap_reflectivity_50deginci_isabel_example

figure(2)
set(gcf,'Position',[0 0 1200 300])
v = -16:2:16;
contourf(yc1,zc1,w2',v)
ylabel('Height (km)');
xlabel('Along Track (km)');
title(['Start:  ',datestr(time(1)),'  End:  ',datestr(time(endd))])
colorbar;

set(gcf,'PaperPositionMode','auto');	%sets printed figure to be same size as window figure.
set(gcf, 'PaperUnits', 'inches');  
set(gcf, 'PaperSize', [11 8.5]);   

print('-dpdf',[datestr(time(1)) 'vertvel-snr'])
%close(2)
%print -dpng iwrap_vertical_velocity_50deginci_isabel_example
return;



%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
%*-0-.nc files, outer beam
%tau = 50.00     %degrees
%*-1-.nc files
%tau = 44.10     %degrees
%*-2-.nc files
%tau = 39.10     %degrees
%*-3-.nc files, inner beam
%tau = 34.50     %degrees




%compute x,y sfc position (relative to aircraft) of each binned measurement, meters
%radius = alt.*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

%for i=1:399                                       
%x = reshape(ilon(:,i,:),32,108);                          
%y = reshape(ilat(:,i,:),32,108);                          
%scatter(x(:,108),y(:,108));
%hold on 
%end
%plot((alon(1,:)+360),alat(1,:),'-k','linewidth',4)  





return;

%m = 0;
%for i=1:32
%    for j=1:399
%        for k=1:115
%            if (xr(i,j,k) >= -0.1 & xr(i,j,k) <= 0.1 & yr(i,j,k) >= 9.9 & yr(i,j,k) <= 10.1 & zr(i,j,k) >= 0.585 & zr(i,j,k) <= 0.615)
%                m = m + 1
%            end;
%        end
%    end
%end

%return;

ix = 0;
for x = -3:0.2:3
    ix = ix + 1
    iy = 0;
    for y = -2:0.2:45
        iy = iy + 1;
        iz = 0;
        for z = 0.030:0.030:2.40
            iz = iz + 1;
   

%Masking out data in grid cell of interest
mask = xr > (x-0.1) & xr < (x+0.1) & yr > (y-0.1) & yr < (y+0.1) & zr > (z-0.015) & zr < (z+0.015);

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

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

%for i=1:32
%    x3 = reshape(xr(i,:,:),size(xr,2)*size(xr,3),1);
%    y3 = reshape(yr(i,:,:),size(xr,2)*size(xr,3),1);
%    z3 = reshape(zr(i,:,:),size(xr,2)*size(xr,3),1);
%    scatter3(x3,y3,z3);
%    hold on
%end

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

return;

%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]);   


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
return;






%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);


return;





ires = 0.10/60;	%initial pulse resolution (will be multiplied by 60 later), km
x = [-1:ires:1];y = x;		%Fine resolution
[x,y]=meshgrid(x,y);
amp = 100	%initial amplitude factor
wide = -550	%initial width factor
wide = -50	%initial width factor

for j=1:10	%width loop
z5 = amp*exp(wide*j*2.25*(x.^2 + y.^2));
cen = (length(x)-1)/2;
cen = cen + 1;
x(cen,cen)	%verifies center
efold = amp/exp(1);	%e-folding amplitude
efold = 0.1;

diam(j) = length(find(z5(cen,:)>=efold))*ires*60;	%diameter of heating pulse
%for i=10:10	%resolution loop
i=10;

z1 = amp*exp(wide*j*2.25*((x-0.416667).^2 + (y-0.00).^2));

%z1 = 100*exp(-10*((x-0.40).^2 + (y-0.40).^2));
%z1 = 100*exp(-150*((x-0.00).^2 + (y-0.00).^2));

deg(i) = i*3		%resolution degradation factor
x2 = [-1:ires*deg(i):1];y2 = x2;	%Coarse resolution
[x2,y2]=meshgrid(x2,y2);
z2 = amp*exp(wide*j*2.25*((x2-0.416667).^2 + (y2-0.00).^2));

%0.5 percent of max heating
%m1 = max(max(z1))*0.005	
%m2 = max(max(z2))*0.005	
%pts1 = length(find(z1>m1))/(length(z1).^2);	%fraction of grid covered by fine resolution bubble
%pts2 = length(find(z2>m2))/(length(z2).^2);	%fraction of grid covered by coarse resolution bubble

%ratio(i,j) = (length(z1).^2)/(length(z2).^2);


%Playing around with plotting...
%figure(1);contourf(x*60,y*60,z1)
%print -dpng test
%return;
%plot(x,z1);return;
%figure(2);contourf(x2*60,y2*60,z2);return;
%hold on;plot(x2,z2,'r');xlim([-0.5 0.5]);return;
%[x,y]=meshgrid(x,y);
%z1 = 100*exp(-50*((x-0.40).^2 + (y-0.40).^2));
%[x,y]=meshgrid(x,y);
%x2 = [-1:0.20:1];y2 = x2;
%[x2,y2]=meshgrid(x2,y2);
%z2 = 100*exp(-50*((x2-0.40).^2 + (y2-0.40).^2));
%contourf(x,y,z1);figure(2);contourf(x2,y2,z2);
%return;
%z2 = 50*exp(-50*((x-0.40).^2 + (y-0.40).^2));
%view(-167,65);zoom(1.125)
%Below fixes the colorbar range...
%clim([0 100]);colorbar

figure(1);surf(x*60,y*60,z1,'edgecolor','none');shading interp
view(-20,60);axis([-60 60 -60 60 0 100]);colorbar;caxis([0 100]);
print -dpng surf1

figure(2);surf(x2*60,y2*60,z2,'edgecolor','none');shading interp;
view(-20,60);axis([-60 60 -60 60 0 100]);colorbar;caxis([0 100]);
print -dpng surf2
%figure(2);contourf(z);hold on;contourf(z2);

return;


%Open binary file for reading ('r') with the precision (i.e. little-endian) of the current system ('n')
%fid = fopen('lh_composite_dropsondes.dat','r','n');
%test = fread(fid,'float');
%new = reshape(test,60,60,20,10,3);
%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);

%x = [-59:2:59];y=x;
%[x,y]=meshgrid(x,y);
%Generate 2-D cylindrical grid
%[r,t] = meshgrid(2:2:60,pi/36:pi/36:2*pi);
%[nx,ny]=pol2cart(t,r);
%zi = interp2(x,y,lh(:,:,5,1),nx,ny,'linear');
%mean1 = nanmean(zi,1)

%cx = [-59:4:59];cy=cx;
%[cx,cy]=meshgrid(cx,cy);
%zic = interp2(x,y,lh(:,:,5,1),cx,cy,'linear');
%contourf(cx,cy,zic.')

%Generate 2-D cylindrical grid
%zi2 = interp2(cx,cy,zic,nx,ny,'linear');
%mean2 = nanmean(zi2,1)

%cx = [-59:8:59];cy=cx;
%[cx,cy]=meshgrid(cx,cy);
%zic = interp2(x,y,lh(:,:,5,1),cx,cy,'linear');

%zi3 = interp2(cx,cy,zic,nx,ny,'linear');
%mean3 = nanmean(zi3,1)

%cx = [-59:10:59];cy=cx;
%[cx,cy]=meshgrid(cx,cy);
%zic = interp2(x,y,lh(:,:,5,1),cx,cy,'linear');

%zi4 = interp2(cx,cy,zic,nx,ny,'linear');
%mean4 = nanmean(zi4,1)



%plot([2:2:60],mean1);hold on;plot([2:2:60],mean2,'r');plot([2:2:60],mean3,'g');plot([2:2:60],mean4,'k')

%return;

%Generate 2-D cylindrical grid
%[r,t] = meshgrid(2:2:60,pi/36:pi/36:2*pi);
[r,t] = meshgrid(0.50:0.50:60,pi/144:pi/144:2*pi);
[nx,ny]=pol2cart(t,r);

zi = interp2(x*60,y*60,z1,nx,ny,'linear');
mean1 = nanmean(zi,1);

zi2 = interp2(x2*60,y2*60,z2,nx,ny,'linear');
mean2 = nanmean(zi2,1);

per = ((max(mean1) - max(mean2))/max(mean1))*100;

%plot([0.5:0.5:60],mean1);hold on;plot([0.5:0.5:60],mean2,'r');
%max(mean1)
%max(mean2)
sp(i,j) = per;
%end
end

contourf(deg,diam,sp);
%ylim([5 40]);
%caxis([-100 100]);colorbar;
colorbar
xlabel('Degradation Factor');
ylabel('Diameter of Heat Pulse (km)');
title('Initial Resolution = 100 m');

%length(find(sp>0));
%figure(2);plot(2:2:60,mean1);hold on;plot(2:2:60,mean2,'r');
%contourf(nx,ny,zi)
%[nr,nt]=cart2pol(nx,ny);
%z = rand(73,30);
