%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;
a=1;
%for a=1:20
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, 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)]);
%time = datestr(test);


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


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

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

%missing = Ze == 0.0;
Ze(missing) = NaN;     

dbz = 10*log10(Ze);
dbz = dbz(ii:ff,:);

%%dbz = shiftdim(dbz(:,ii:ff),1);
dbz = reshape(dbz,32,pp/32,183);

return;

avgf = fspecial('average');
pass = squeeze(dbz(1,:,:));
new = imfilter(pass,avgf);

return;

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

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

fix1 = bsxfun(@minus,vel,vcorr);

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:183]*30)/1000;     %km
roll = nav(:,:,11);
pitch = nav(:,:,10);
drift = nav(:,:,8);
tau = 50.00;

maptrack;

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

return;

%Reflectivity processing...
missing = Ze <= 0.0;      %if inequality is true, missing is 1; false, missing is 0.
Ze(missing) = NaN;     
dbz = 10*log10(Ze);

dbz = shiftdim(dbz(:,ii:ff),1);
dbz = smooth3(reshape(dbz,32,599,183));
%dbz = rot90(dbz(:,10:100));
v = 30:2:60;
return;
%to get correct plot transpose array (using "'") or rotate 2D matrix by 90 degrees
%Also need to reverse y axis direction
contourf(dbz',v)
set(gca,'YDir','reverse')
%dbz = rot90(dbz(:,10:100));
%End Reflectivity processing

return;

%Doppler Velocity processing...
vel = shiftdim(vel(:,ii:ff),1);

thetap = thetap(ii:ff);
thetap = reshape(thetap,32,599);

velz = bsxfun(@times,vel,cos(thetap*d2r));

return;

vcorr = vcorr(ii:ff);
vcorr = reshape(vcorr,32,599);

fix = bsxfun(@minus,vel,vcorr);

svel = svel(ii:ff);
svel = reshape(svel,32,599);






return;


%Note:  variables are times (seconds) relative to 1-jan-0000 00:00:00
sec1970 = datenum('1-jan-1970 00:00:00')*24*3600;   
start = datenum('12-sep-2003 19:01:14')*24*3600;
endd = datenum('12-sep-2003 19:07:53')*24*3600;

%Variables below are unix time stamps
start = start - sec1970;
endd = endd - sec1970;

%Selecting period of time that is consistent with other files...
mask = unixt >= start & unixt <= endd;


return;

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

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