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_geometry.txt');
g=textscan(fid,'%s');

fid = fopen('passes_cband_reflectivity.txt');
z=textscan(fid,'%s');

fid = fopen('passes_cband_doppler.txt');
k=textscan(fid,'%s');

ncid = netcdf.open('nav030912H.nc','NC_NOWRITE');

%Variable indices...
%0 = time, 11 = ground speed of aircraft

atime = netcdf.getVar(ncid,0);
gs = netcdf.getVar(ncid,11);

%for i=1:size(gs)
%atime2 = datenum([1970 1 1 0 0 atime(i)]);
%atime2 = datestr(atime2);
%end

%return;

cnt = 0;
a=1;
%for a=1:20
path = fullfile(g{1}{a})
ncid = netcdf.open(path,'NC_NOWRITE');

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


return;

igs = interp1(atime,gs,unixt,'linear');     %ground speed of aircraft interpolated to scan pattern

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);
alt = 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;

%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


%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 = alt./cos(iang*d2r);
sfc_rg = round(mean(mean(rsfc))/30);		%mean surface range gate	


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


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


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:398
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
%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

Taut = asin(bsxfun(@rdivide,Yt,reshape(range,1,1,size(range,2))));

%Get reflectivity and Doppler data

path = fullfile(k{1}{a})
ncid = netcdf.open(path,'NC_NOWRITE');

%Variable indices...Note that names should be swapped in the .nc files
%0 = vertical component of doppler velocity
%1 = horizontal component of doppler velocity
DVz = netcdf.getVar(ncid,0); 	%units m/s
DVh = netcdf.getVar(ncid,1);    %units m/s
%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

iang = reshape(iang,1,size(iang,1),size(iang,2));

DV = bsxfun(@rdivide,DVz,cos(iang*d2r));

niang = iang(1,:,200);
nazi = azi(:,200);nazi = reshape(nazi,1,32);

ndvh = DVh(112,:,200);
ndvz = DVz(112,:,200);

th = ndvh.*sin(niang*d2r).*cos(nazi*d2r) 

%+ ndvz.*cos(niang*d2r);

return;

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

in = bsxfun(@times,igs,sin(Taut)) + DV;

nDV = double(reshape(in,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,nDV);
DVg = F(xc,yc,zc);
DVg = smooth3(DVg);

return;


%path = fullfile('/mercury/guimond/GSFC/IWRAP/20030912/volume_reflectivity/',z{1}{16})
path = fullfile(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);
