function [vel_corrected,x_grid,y_grid,z_grid]=cressman(advection,first_file)
volume_size=7;
%first_file=1;
hurricane='frances';
rep='meanless\vert';
path=['C:\Documents and Settings\sylvie lorsolo\My Documents\research\',hurricane,'\data\modified\',rep,'\'];
%path2=['C:\Documents and Settings\sylvie lorsolo\My Documents\research\',hurricane,'\data\total_wind\'];
VAD_data(600,1)=0;
SO=150;
ROX=150;
ROY=150;
ROZ=150;
max_grid_horiz=14000;
max_grid_vert=2000;
delta_x=50;
delta_y=50;
delta_z=25;
x_grid=-2500:delta_x:15000;
y_grid=-2500:delta_y:15000;
z_grid=0:delta_z:max_grid_vert;
delta_profile=50;
z_grid_profile=0:delta_profile:((floor((max(z_grid)+ROZ)/delta_profile)+1)*delta_profile);
z_advection=300;

[x_3dgrid,y_3dgrid,z_3dgrid]=meshgrid(x_grid,y_grid,z_grid);

[x_size,y_size,z_size]=size(x_3dgrid);
w_tot(x_size,y_size,z_size)=0;
w_by_vel_tot(x_size,y_size,z_size)=0;
vel_corrected(x_size,y_size,z_size)=0;
w_tot_profile(length(z_grid_profile))=0;
count(x_size,y_size,z_size)=0;
w_by_wd_tot(length(z_grid_profile))=0;
w_by_ws_tot(length(z_grid_profile))=0;






%Obtain files
%total_wind_files=dir([path2,'*.dat']);



if advection==1 | advection == 2
    first_file_vad=1;
    last_file=first_file_vad + volume_size-1;
    t0=sweeptime(total_wind_files(first_file_vad).name);
    t_final=sweeptime(total_wind_files(last_file).name);
    delta_t=etime(t_final,t0);
    median_time=floor((delta_t)/2);
    for fil=first_file_vad:last_file
        vad_file=total_wind_files(fil).name;
        VAD_data=load([path2,vad_file]);
        el=str2num(vad_file(30:32)); % works for el < 10 deg
        rge=33.35:66.7:(length(VAD_data)-1)*66.7+33.35;
        z_obs=rge*sin(deg2rad(el));
        
        
        for r=1:length(VAD_data)
            zminp=z_obs(r)-ROZ;
            zmaxp=z_obs(r)+ROZ;
            u_z=find(z_grid_profile >= zminp & z_grid_profile <= zmaxp);
            zminp=min(u_z);
            zmaxp=max(u_z);
            for k=zminp:zmaxp
                d_profile=abs(z_grid_profile(k)-z_obs(r));
                if d_profile <= ROZ & VAD_data(r,1)<100 & VAD_data(r,2)<=360 & VAD_data(r,1)>=0 & VAD_data(r,2)>=0;
                    w_profile=(ROZ^2 - d_profile^2)/(ROZ^2 + d_profile^2);
                    w_tot_profile(k)=w_tot_profile(k)+w_profile;
                    w_by_wd_tot(k)=w_by_wd_tot(k)+w_profile*VAD_data(r,2); %wind direction
                    w_by_ws_tot(k)=w_by_ws_tot(k)+w_profile*VAD_data(r,1);
                end
            end
            
        end
        
        u_0=find(w_tot_profile==0);
        wd_corrected(u_0)=NaN;
        ws_corrected(u_0)=NaN;
        
        u_non0=find(w_tot_profile~=0);
        wd_corrected(u_non0)=w_by_wd_tot(u_non0)./w_tot_profile(u_non0);
        ws_corrected(u_non0)=w_by_ws_tot(u_non0)./w_tot_profile(u_non0);
    end
    
    if advection == 2
        u_z=find(z_grid_profile >= z_advection-ROZ & z_grid_profile <= z_advection + ROZ);
        zminp=min(u_z);
        zmaxp=max(u_z);
        ws_mean_layer=0;
        wd_mean_layer=0;
        for k=zminp:zmaxp
            ws_mean_layer=ws_mean_layer+ws_corrected(k);
            wd_mean_layer=wd_mean_layer+ws_corrected(k);
        end
        ws_mean_layer=ws_mean_layer/length(u_z);
        wd_mean_layer=wd_mean_layer/length(u_z);
    end
end

%--------------------------------------------GRIDDING ALGORITHM--------------------------------%%%%%% 
radar_files=dir([path,'*.dat']); 

last_file=first_file + volume_size-1;

for fil=first_file:last_file
    sweepfile=radar_files(fil).name;
    radical=sweepfile(1:length(sweepfile)-4);
    %VAD_data=load([path2,radical,'_wind.dat']);
    data=load([path,sweepfile]);
    az_N=data(length(data),:);
    vel=data(1:length(VAD_data),:);
    el=str2num(sweepfile(30:32));
    
    
    rge=33.35:66.7:(length(VAD_data)-1)*66.7+33.35;
    z_obs=rge*sin(deg2rad(el));
    
    az=90-az_N;
    az_x=deg2rad(az);
    rge2d=rge'*ones(1,length(az));
    az2d=ones(length(VAD_data),1)*az_x;
    el=el*ones(length(VAD_data),length(az));
    [x,y,z]=sph2cart(az2d,deg2rad(el),rge2d);
    
    for r=1:length(VAD_data)
        if advection == 1
            t_temp=sweeptime(sweepfile);
            sec_diff=etime(t_temp,t0);
            time_diff=median_time-sec_diff;
            u_zp=find(z_grid_profile> (z(r,1)-delta_profile/2) & z_grid_profile<= (z(r,1)+delta_profile/2));
            if  VAD_data(r,1)<100 & VAD_data(r,2)<=360 & VAD_data(r,1)>0 & VAD_data(r,2)>=0 & length(u_zp)~=0
                delta_d=time_diff*ws_corrected(u_zp(1));
                wd_x=90-(wd_corrected(u_zp(1))-180);
                x_adv_corr(r,:)=x(r,:)+ cos(deg2rad(wd_x))*delta_d;
                y_adv_corr(r,:)=y(r,:)+ sin(deg2rad(wd_x))*delta_d;
            else
                x_adv_corr(r,:)=x(r,:);
                y_adv_corr(r,:)=y(r,:);
            end
        elseif advection == 2
            t_temp=sweeptime(sweepfile);
            sec_diff=etime(t_temp,t0);
            time_diff=median_time-sec_diff;
            if  VAD_data(r,1)<100 & VAD_data(r,2)<=360 & VAD_data(r,1)>0 & VAD_data(r,2)>=0
                delta_d=time_diff*ws_mean_layer;
                wd_x=90-(wd_mean_layer-180);
                x_adv_corr(r,:)=x(r,:)+ cos(deg2rad(wd_x))*delta_d;
                y_adv_corr(r,:)=y(r,:)+ sin(deg2rad(wd_x))*delta_d;
            else
                x_adv_corr(r,:)=x(r,:);
                y_adv_corr(r,:)=y(r,:);
            end
        elseif advection == 3
            t_temp=sweeptime(sweepfile);
            sec_diff=etime(t_temp,t0);
            time_diff=median_time-sec_diff;
            if  VAD_data(r,1)<100 & VAD_data(r,2)<=360 & VAD_data(r,1)>0 & VAD_data(r,2)>=0
                ws_mean_layer=7.5;
                wd_mean_layer=31;
                wd_x=90-(wd_mean_layer-180);
                delta_d=time_diff*ws_mean_layer;
                x_adv_corr(r,:)=x(r,:)+ cos(deg2rad(wd_x))*delta_d;
                y_adv_corr(r,:)=y(r,:)+ sin(deg2rad(wd_x))*delta_d;
            else
                x_adv_corr(r,:)=x(r,:);
                y_adv_corr(r,:)=y(r,:);
            end
        else
            x_adv_corr(r,:)=x(r,:);
            y_adv_corr(r,:)=y(r,:);
        end
        
        for azi=1:length(az)
            xmin=x_adv_corr(r,azi)-ROX;   
            xmax=x_adv_corr(r,azi)+ROX;
            ymin=y_adv_corr(r,azi)-ROY;
            ymax=y_adv_corr(r,azi)+ROY;
            zmin=z(r,azi)-ROZ;
            zmax=z(r,azi)+ROZ;
            xfind=find(x_grid>=xmin & x_grid<=xmax);
            yfind=find(y_grid>=ymin & y_grid<=ymax);
            zfind=find(z_grid>=zmin & z_grid<=zmax);
            xmin=min(xfind);
            xmax=max(xfind);
            ymin=min(yfind);
            ymax=max(yfind);
            zmin=min(zfind);
            zmax=max(zfind);
            
            for u=xmin:xmax
                for v=ymin:ymax
                    for ww=zmin:zmax                     
                        dist=((y_3dgrid(v,u,ww)-y_adv_corr(r,azi))^2+(x_3dgrid(v,u,ww)-x_adv_corr(r,azi))^2+(z_3dgrid(v,u,ww)-z(r,azi))^2)^0.5;
                        if dist<=SO
                            w=(SO^2-dist^2)/(SO^2+dist^2);
                            if vel(r,azi)<100 & vel(r,azi)>-100 & isnan(vel(r,azi))~=1
                                w_tot(v,u,ww)=w_tot(v,u,ww)+w;
                                w_by_vel_tot(v,u,ww)=w_by_vel_tot(v,u,ww)+w*vel(r,azi);
                                count(v,u,ww)=count(v,u,ww)+1;
                            end
                        end                            
                    end
                end
            end  
        end
    end
end
for u=1:x_size %iterate through final cartesian grid and complete Creesman
    for v=1:y_size
        for ww=1:z_size
            if w_tot(v,u,ww)==0
                vel_corrected(v,u,ww)=NaN;
            else vel_corrected(v,u,ww)=w_by_vel_tot(v,u,ww)/w_tot(v,u,ww);
            end
        end
    end
end   

u_threshold=find(count<20);
%vel_corrected(u_threshold)=NaN;
% figure
% pcolor(x_grid,y_grid,vel_corrected(:,:,5))
% shading interp
% hold on
% h7=pcolor(x_grid,y_grid,vel_corrected(:,:,7));
% shading interp
% zd7=get(h7,'ZData');
% zd7=2+zd7;
% set(h7,'zdata',zd7)
% h10=pcolor(x_grid,y_grid,vel_corrected(:,:,10));
% zd10=get(h10,'ZData');
% zd10=5+zd10;
% set(h10,'zdata',zd10)
% shading interp
%[x_3dgrid,y_3dgrid,z_3dgrid]=meshgrid(xc,yc,zc);
% x1=-10000:50:10000;
% y1=-10000:50:10000;
% hsp = surf(x1,y1,zeros(401));
% rotate(hsp,[1,0,0],90)
% rotate(hsp,[0,0,1],135)
% xd = get(hsp,'XData');
% yd = get(hsp,'YData');
% zd = get(hsp,'ZData');
% xd1=xd+5000;
% yd1=yd+5000;
% h=slice(x_3dgrid,y_3dgrid,z_3dgrid,vel,xd1,yd1,zd);
% x_section = get(h,'CData');
% pcolor(xd1(1,:),zd(:,1),x_section)
% shading interp
%to display only the right upper corner, for x and y positive:
% u=find(xc(1,:)>=0);
% xc2=xc(u);
% yc2=yc(u);
% vel2=vel(u,u,:);
% pcolor(xc2,yc2,vel2(:,:,10))

vel=vel_corrected;
xc=x_grid;
yc=y_grid;
zc=z_grid;
sweepfile=radar_files(first_file).name;
time=sweepfile(12:17);
save(['C:\Documents and Settings\sylvie lorsolo\My Documents\research\frances\data\modified\meanless\vert\gridded_data\vol',time],'vel','xc','yc','zc');