function [x,y,z,u_obs,v_obs,w_obs, p_position_x,p_position_y]=location_3d( u_grid, v_grid, w_grid, inc_iwrap_t, az_N, Vp, trackp, alt_p, T, delta_d_plane)
%all input should be (N,1) except the wind component and the cartesian grid
tic

% if a predefined wind field and cartesian grid is entered predef=1
% ortherwise it is 0 and it will be defined in the program
predef=1;


if nargin < 3
    
    T=20;
    time_rotation=1;
    delta_az=4;
    delta_t= (delta_az*time_rotation)/360;
    N=T/delta_t
    Vp=120*ones(N,1);
    delta_d_plane=delta_t*Vp;
    az_N=0:delta_az:356;
    az_a=90-az_N;
    u_neg=find(az_a<0);
    az_a(u_neg)=az_a(u_neg)+360;
    az=(az_a)'*ones(1,(T/time_rotation));
    az=reshape(az,length(az)*T/time_rotation,1);
    inc_iwrap_t=29.49*ones(N,1);
    trackp=60*ones(N,1);
    alt_p=1500*ones(N,1);;
else
    az=90-az_N;
    u_neg=find(az<0);
    az(u_neg)=az(u_neg)+360;
    N=length(az_N);
end

rge=15:30:56*30+15;
HR=rge'*sin(deg2rad(inc_iwrap_t'));
p_position_x(1,1)=0;
p_position_y(1,1)=600;

for i=1:N-1
    p_position_x(i+1,1)=delta_d_plane(i+1).*sin(deg2rad(trackp(i+1))) + p_position_x(i,1);
    p_position_y(i+1,1)=delta_d_plane(i+1).*cos(deg2rad(trackp(i+1))) + p_position_y(i,1);  
end

if predef==0
delta_x=150;
delta_y=150;
delta_z=30;
x_grid=0:delta_x:1000;  
x_grid=[-1000:delta_x:-150 x_grid];
y_grid=x_grid;
z_grid=0:50:1500;
[x_3dgrid,y_3dgrid,z_3dgrid]=meshgrid(x_grid,y_grid,z_grid);
else
    [x_3dgrid,y_3dgrid,z_3dgrid,u_grid,v_grid,w_grid]=roll_windfield;
end





%[th,r]=meshgrid(deg2rad(az'),HR);
[x,y]=pol2cart(ones(length(rge),1)*deg2rad(az)',HR);
x=x+ones(length(rge),1)*p_position_x';
y=y+ones(length(rge),1)*p_position_y';% it should be displacement if the fist plane position is not at 0
z=ones(length(rge),1)*alt_p'-rge'*cos(deg2rad(inc_iwrap_t'));

u_obs=zeros(size(x));
v_obs=zeros(size(x));
w_obs=zeros(size(x));

% u_obs=interp3(x_grid,y_grid,z_grid,u_grid,x,y,z);
% v_obs=interp3(x_grid,y_grid,z_grid,v_grid,x,y,z);
% w_obs=interp3(x_grid,y_grid,z_grid,w_grid,x,y,z);
disp('about to start cress')
[x_size,y_size]=size(x);
[d,u_obs,v_obs,w_obs]=cres1(x,y,z, x_3dgrid,y_3dgrid,z_3dgrid, u_grid, v_grid, w_grid);
u999=find(u_obs==-999.0);
u_obs(u999)=NaN;
v_obs(u999)=NaN;
w_obs(u999)=NaN;
% for i=1:x_size
%     for j=1:5%:y_size
%         d=sqrt((x_3dgrid-x(i,j)).^2 + (y_3dgrid-y(i,j)).^2 + (z_3dgrid-z(i,j)).^2 );
%         ud=find(d<=RO);
%         if length(ud)>2
%             weight=(RO^2 - d(ud).^2)./(RO^2 + d(ud).^2);
%             u_obs(i,j)=sum(weight.*u_grid(ud))/sum(weight);
%             v_obs(i,j)=sum(weight.*v_grid(ud))/sum(weight);
%             w_obs(i,j)=sum(weight.*w_grid(ud))/sum(weight);
%         else
%             u_obs(i,j)=NaN;
%             v_obs(i,j)=NaN;
%             w_obs(i,j)=NaN;
%             
%             
%         end
%     end
% end

toc

