test = load('/mimosa/jbeck/home/RIP4_v4.3/test.out');

for i=1:(length(test)/181)
    trajs(i,:,:)=test((((i-1)*181)+1):(i*181),:);
end

%Hours 1 to 3 in decimal numbers
times=1:(1/60):3;

%Find the location of the times corresponding to the trajectories with data

for i=1:length(times)
    for j=1:length(trajs(1,:,1))
        if abs(trajs(1,j,1)-times(i)) < 0.0005
            trajsdata(:,j,:)=trajs(:,j,:);
        end
    end
end

nonzero = find(trajsdata(1,:,1)~=0);

for i=1:length(nonzero)
    trdata(:,i,:) = trajsdata(:,nonzero(i),:);
end

time = trdata(1,:,1)*60;
time = round(time);
time = time/60;
wrftime = time2str(time);

files=dir('/media/wrfdata/control/workspaces/wrfout*');

%Barnes interpolation scheme to compute data values along trajectories by
%interpolating from model grid points values

%Sphere of influence in meters
SO=500;

%Radius of influence along each axis
ROX=SO;
ROY=SO;
ROZ=SO;

%Horizontal domain in grid points
max_grid_horiz=50000;

%Grid spacing
delta_x=250;
delta_y=250;

%Define the grid
x_grid=0:delta_x:max_grid_horiz;
y_grid=0:delta_y:max_grid_horiz;
z_grid=[100;300;500;700;950;1200;1400;1600;1800;2000;2300;2500;2700;3000;3200;3400;3700;3900;4200;4400;4700;4900;5200;5500;5700;6000;6300;6600;6900;7200];

%Create the grid
[x_3dgrid,y_3dgrid,z_3dgrid]=meshgrid(x_grid,y_grid,z_grid);

%Define the size of the grid
[x_size,y_size,z_size]=size(x_3dgrid);

%Initialize the barnes weighting variables to 0
w_tot(length(trdata(:,1,1)),length(trdata(1,:,1)))=0;
w_by_var_tot(length(trdata(:,1,1)),length(trdata(1,:,1)))=0;
var_corrected(length(trdata(:,1,1)),length(trdata(1,:,1)))=0;

%--------------------------------------------GRIDDING ALGORITHM--------------------------------%%%%%% 
for e=1:length(x_size)

   for f=1:length(y_size)
       
       for g=1:length(z_size)
       
       %Load data files
        
       load(['/media/wrfdata/control/workspaces/wrfout_d01_0001-01-01_',wrftime(f,:),'-box.mat']);
       
       
                %xmin=(trajs(e,nonzero(f),2)*250)-ROX;
                xmin=x_grid(e)-ROX;
                xmax=x_grid(e)-ROX;
                ymin=y_grid(f)-ROY;
                ymax=y_grid(f)-ROY;
                zmin=z_grid(g)-ROZ;
                zmax=z_grid(g)-ROZ;
                xfind=find(trdata(1,:,2)*250>=xmin & trdata(1,:,2)*250<=xmax);
                yfind=find(trdata(1,:,3)*250>=ymin & trdata(1,:,3)*250<=ymax);
                zfind=find(trdata(1,:,4)*1000>=zmin & trdata(1,:,4)*1000<=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(e,f,g)-(trdata(1,v,3)*250))^2+(x_3dgrid(e,f,g)-(trdata(1,u,2)*250))^2+(z_3dgrid(e,f,g)-(trdata(1,ww,4)*1000))^2)^0.5;
                            if dist<=SO
                                %w=(SO^2-dist^2)/(SO^2+dist^2);%Cressman
                                w=exp(-(((dist)/(SO))^2)); %Barnes
                                w_tot(u,v,ww)=w_tot(u,v,ww)+w;
                                w_by_var_tot(u,v,ww)=w_by_var_tot(u,v,ww)+w*VERTVORTfinal(ww,v,u);                                  
                            end
                        end
                    end
                end
                
                for u=size(length(trajs(:,1,1)))
                
                    for v=size(length(trajs(:,1,1)
                if w_tot(u,v,ww)==0 
                    var_corrected(u,v,ww)=NaN;
                else var_corrected(u,v,ww)=w_by_var_tot(u,v,ww)/w_tot(u,v,ww);
                end
                
      end
   end
end