%small experiment along the Lebanese coast

%floaters along east coast
close all;
clear all; 
write_ini=1;


Nlat=58; %Phi
Nlong=87; %Lambda
load 'meshgrid2.dat'; 

GR_long=reshape(meshgrid2(:,1),Nlong,Nlat); %will give a Nlong by Nlat matrix 
%containing the long of the points 
GR_lat=reshape(meshgrid2(:,2),Nlong,Nlat); %will give a Nlong by Nlat matrix 
%containing the lat of the points 

del=GR_long(2)-GR_long(1);


%for plotting
load ../obs_float/uzero.dat;
load ../obs_float/vzero.dat;
uu=reshape(uzero(:,4),58,87);
vv=reshape(vzero(:,4),58,87);


%5 cities in lebanon: Tripoli, Beirut, Byblos, Tyre and Nacoura


%translate from long and lat in degrees to mesh points
[cc1 v1]=find(abs(GR_long-35.53)<del/2); 
[u2 cc2]=find(abs(GR_lat-34.41)<del/2); 
f1_1=cc1(1)+((35.53-GR_long(cc1(1)))./del);
f1_2=cc2(1)+((34.41-GR_lat(1,cc2(1)))./del);


[bb1 v1]=find(abs(GR_long-35.38)<del/2); 
[u2 bb2]=find(abs(GR_lat-34.12)<del/2); 
f2_1=bb1(1)+((35.38-GR_long(bb1(1)))./del);
f2_2=bb2(1)+((34.12-GR_lat(1,bb2(1)))./del);


[bb1 v1]=find(abs(GR_long-35.253)<del/2); 
[u2 bb2]=find(abs(GR_lat-33.88)<del/2); 
f3_1=bb1(1)+((35.253-GR_long(bb1(1)))./del);
f3_2=bb2(1)+((33.88-GR_lat(1,bb2(1)))./del);


[bb1 v1]=find(abs(GR_long-35.18)<del/2); 
[u2 bb2]=find(abs(GR_lat-33.6)<del/2); 
f4_1=bb1(1)+((35.18-GR_long(bb1(1)))./del);
f4_2=bb2(1)+((33.6-GR_lat(1,bb2(1)))./del);

[bb1 v1]=find(abs(GR_long-34.98)<del/2); 
[u2 bb2]=find(abs(GR_lat-33.27)<del/2); 
f5_1=bb1(1)+((34.98-GR_long(bb1(1)))./del);
f5_2=bb2(1)+((33.27-GR_lat(1,bb2(1)))./del);


%store in an array
N_f=5;%number of drifters
FF=zeros(N_f,2);
FF(1,:)=[f1_1 f1_2];FF(2,:)=[f2_1 f2_2];FF(3,:)=[f3_1 f3_2];FF(4,:)=[f4_1 f4_2];
FF(5,:)=[f5_1 f5_2];



if (write_ini==1)
    FF=FF-1; %yao to matlab conversion (i=0,1, in YAO)
    
    
    %Put in YAO format
    
    xy_vec=ones(N_f*2,1); xy_vec(2:2:end)=2;
    time_vec=ones(N_f*2,1);
    
    %floater id
    Nf_vec=[1; 1];
    for j=1:N_f-1;
        Nf_vec=[Nf_vec; j+1 ;j+1];
    end
    
    %the positions
    pos_vec=[FF(1,1); FF(1,2)];
    for j=2:N_f;
        pos_vec=[pos_vec; FF(j,1); FF(j,2)];
    end
    
    
    
    Yao_F=[ floor(xy_vec)  floor(time_vec) floor(Nf_vec) pos_vec];
    
    
    dlmwrite('../obs_float/init5.dat',Yao_F,'\t');
    
end;


figure;
quiver(GR_long(45:end,1).',GR_lat(1,:),uu(:,45:end),vv(:,45:end))
hold on;

for j=1:N_f

F11=floor(FF(j,1));F22=floor(FF(j,2));
alphaLong=(FF(j,1)-F11);
alphaLat=(FF(j,2)-F22);
plot(   GR_long(F11,F22)+alphaLong*del, GR_lat(F11,F22)+alphaLat*del,'-xr');

end

xlabel('Longitudes en degrees'); ylabel('Latitudes en degrees'); 
