% Task 5.0: Ditch and Channel filling 

% 7 sites of interest and three transects at each = 21 transects total 
% Present: 1 vegetation scheme x 4 storms x 2 tide conditions = 8 results per transect 
% In 2100: 1 vegetation scheme x 4 storms x 2 tide conditions x 1 SLR scenario  = 8 results per transect
% -	16 results per transect for all transects 
% -	Mean vegetation conditions will be used. 
% -	The storms are the four already selected. 
% -	The tide conditions correspond to low and high tide. 
% -	The SLR scenario is high (RCP 8.5) projected to the year 2100.
tic
area1=[5,23,30];
area2=[35,14,22];
area3=[31,35,36];
area4= [0,1,2]; %[27,28,29];
area5=[30,14,32];
area6=[29,38,21];
area7=[19,10,6];
AreaList=[area1; area2; area3; area4; area5; area6; area7];
tide_conditions = {'L','H'};
R_WA = nan(1000,21);
R_Elev = nan(1000,21);
R_Dep = nan(1000,21);
R_CDa = nan(1,21);
wave_threshold = 0.05;

for n = 4:7%storm number: 4,5,6,7
    for td = 1:2 
        tide = tide_conditions{td}; %tide condition: 'H' or 'L'
        storm_name = ['C0' num2str(n) '_' tide '_Base']
        output_wave_name = ['WaveAtten_5_3_0_' storm_name];
        output_elv_name = ['TrElev_5_3_0_' storm_name];
        output_dep_name = ['TrDep_5_3_0_' storm_name];
        output_CD_name = ['TrCD_5_3_0_' storm_name];
        
        %Get the wave and depth data from this storm run
        cd(['A:\Newbury_new_runs\New Storm Run Corrected\' storm_name]);
        fidw=qpfopen(['wavm-' storm_name '.dat']);
        H = qpread(fidw,'hsig wave height','griddata',0);
        T = qpread(fidw,'peak wave period','griddata',0);
        %h = qpread(fidw,'water depth','griddata',0);       
        ub = qpread(fidw,'orbital velocity near bottom','griddata',0);      
        lm = qpread(fidw,'mean wave length','griddata',0);    
        
        
        fidh = qpfopen(['trim-' storm_name '.dat']);
        h = qpread(fidh,'water depth','griddata',0);      
        bl = qpread(fidh,'bed level in water level points','griddata',0);
                  
        cd('C:\Users\mrfoster\OneDrive - University of New Orleans\Documents\Newbury_Modeling\Transect_points')
        j = 1;
        for An =[1,2,3,4,5,6,7] %1:7 %areas of interest
            disp(['On Area ' num2str(An) ' of Interest'])
            for i = 1:3
                disp(['On Transect ' num2str(i)])
                A = shaperead(['Trans_pt_area' num2str(An) 's.shp'],'Selector',{@(v1) (v1 == AreaList(An,i)), 'ORIG_FID'});
                if An == 3 || An == 4
                    X = A(1).POINT_X;
                    Y = A(1).POINT_Y;
                else
                    X = A(end).POINT_X;
                    Y = A(end).POINT_Y;
                end
                
                dist_2pt = sqrt(((X-H.X).^2)+((Y-H.Y).^2));
                [r,c] = find(dist_2pt == min(min(dist_2pt)));
                mi = find(max(H.Val(1:3,r,c)));
                if isempty(mi)
                    mi = 1;
                end
                H0 = H.Val(mi,r,c);
                ub0 = ub.Val(mi,r,c);
                k0 = (2*pi)./lm.Val(mi,r,c);
                T0 = T.Val(mi,r,c);
                
                dist_2pt = sqrt(((X-h.X).^2)+((Y-h.Y).^2));
                [r,c] = find(dist_2pt == min(min(dist_2pt)));
                h0 = h.Val(mi,r,c); 
                bl0 = bl.Val(r,c);
                
       
                elv = extractfield(A, 'lidar_elev')';
                elv_old = elv;
          

                
              [elv] = fillditches(elv);
              
                
                N = 3093;
                bv = 0.0027;
                hv = 0.27;
                
                [WH,El,Dp,CDa] = transect_wave_attenuation(An, N, bv, ub0, T0, H0, h0, bl0, hv, k0, elv, wave_threshold);
                
                R_WA(1:length(WH),j) = WH;
                R_Elev(1:length(El),j) = El;
                R_Dep(1:length(Dp),j) = Dp;
                R_CDa(j) = CDa;
                j = j+1;
                
            end
        end
        
        cd('A:\Newbury_new_runs\Wave_attenuation_model_results')
        save(output_wave_name,'R_WA') ;
        save(output_elv_name,'R_Elev') ;
        save(output_dep_name,'R_Dep') ;
        save(output_CD_name, 'R_CDa');
        
        
        
    end
end
toc