disp('Indicate if they are data are from Cortadura (area=1) or from Camposoto (area=2)')
area=input('area=');

switch area
    case 1
        site = 'cort';
        load C:\MICORE\Cortadura\cort_allgrids
        load C:\MICORE\Cortadura\clavosXYZ
    case 2
        site = 'camp';
        load C:\MICORE\Camposoto\camp_allgrids
        load C:\MICORE\Camposoto\clavosXYZ
end

if area==1
     disp('Indicate if they are data are from City (area=1) or from Dunes (area=2)')
     area1=input('area=');
     if area1 ==1
         BASE = BASE_CITY;
         REF  = REF_CITY;
     elseif area1==2
         BASE = BASE_DUNES;
         REF  = REF_DUNES;
     end
end





disp('Choose start and end point of profile (c=1) Profile perpenticular to the beach from BASE STATION  (c=2)')
c=input('c=');

switch c
    case 1
        figure
        contour(x,y,ZALL(:,:,1),20)
        hold on 
        scatter(BASE(1),BASE(2),'+k')

        annotation('arrow')
        pause 
        disp('Click origin of the profile x=0,y=0 (landward point)')
        disp('Then click the end of the profile (seaward point)')
        [xi,yi] = ginput(2);
        % Profiles limits
        xori=xi(1);yori=yi(1);
        xend=xi(2);yend=yi(2);
        % Get profiles angle
        alfa=atan2(yend-yori,xend-xori);
        
    case 2
        if area==1
            alfa = -2.5604;
        elseif area ==2
            alfa = -2.7489;
        end
        xori=BASE(1);yori=BASE(2);
        disp('Give the length of the profile dist(m)')
        dist=input('dist=');
        xend=xori+dist*cos(alfa); yend=yori+dist*sin(alfa);;
end

%Calculate points on the profile

disp('Indicate the step of the profile ds(m)')
ds  = input('ds=');

XI  = [xori:ds*cos(alfa):xend];
YI  = [yori:ds*sin(alfa):yend];
lng = sqrt((xori-xend).^2 + (yori-yend).^2);
disp('Indicate how many layer (profiles) to extract)')
k =input('k=');

for i=1:k
    
    ZI(i,:) = interp2(x,y,ZALL(:,:,i),XI,YI);
    
end

% Make Plot
ZI = [[0:ds:lng]' ZI'];
figure
plot(ZI(:,1),ZI(:,2:end));
xlabel('Length (m)')
ylabel('Elevation (m)')




