volume_size=14;
 first_file=1;
hurricane='fran';
rep='meanless\wave1';
path=['C:\Documents and Settings\sylvie lorsolo\My Documents\research\',hurricane,'\data\modified\',rep,'\'];
VAD_data(50,1)=0;
SO=300;
ROX=SO;
ROY=SO;
ROZ=SO;
% max_grid_horiz=14000;
max_grid_vert=2500;
delta_x=100;
delta_y=100;
delta_z=50;
x_grid=-12500:delta_x:12500;                    %for frances_vert_0Z:-9000:delta_x:8000;   ????%-12500:delta_x:2000;    %for frances: -2500:delta_x:15000;         %for isabel_vert:-14500:delta_x:3000;
y_grid=-12500:delta_y:12500;                                        %-9000:delta_x:8000;   ????%-12500:delta_y:2000;                 %-2500:delta_y:15000;                         %-14500:delta_y:3000;
z_grid=0:delta_z:max_grid_vert;



[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;
count(x_size,y_size,z_size)=0;






%--------------------------------------------GRIDDING ALGORITHM--------------------------------%%%%%% 
radar_files=dir([path,'*.dat']); 
last_file=first_file + volume_size-1;
%gamma=1;
%for iter=1:4
%     w_tot(x_size,y_size,z_size)=0;
%     w_by_vel_tot(x_size,y_size,z_size)=0;
%     SO=gamma*SO
%     ROX=SO;
%     ROY=SO;
%     ROZ=SO;
    for fil=first_file:last_file
        sweepfile=radar_files(fil).name;
        data=load([path,sweepfile]);
        [li,col]=size(data);
        az_N=data(li,1:359);
        vel=data(1:li,1:359);
        u_ze=find(vel==0);
        vel(u_ze)=NaN;
        el=str2num(sweepfile(21:23));
        
        
        rge=125:250:49*250+125;%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)
            
            x_adv_corr(r,:)=x(r,:);
            y_adv_corr(r,:)=y(r,:);
            
            
            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);%Cressman
                                %w=exp(-(dist^2)/(2*SO^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
    %gamma=0.8;
    %end
%u_threshold=find(count<15);
%vel_corrected(u_threshold)=NaN;
vel=vel_corrected;
xc=x_grid;
yc=y_grid;
zc=z_grid;
% h10=pcolor(xc,yc,vel(:,:,10));
% shading interp

sweepfile=radar_files(first_file).name;
time=sweepfile(12:17);
save(['C:\Documents and Settings\sylvie lorsolo\My Documents\research\',hurricane,'\data\modified\',rep,'\gridded_data\vol',time],'vel','xc','yc','zc','count');