% Task 5.2: Coverting to tidal flats

% 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
%Three transects were selected out of all the ones generated. These are
%their ORIG_FID: 
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;
task = '2'; 
disp(['OnTask 5 ' num2str(task)])
for v = 4% 1:3 %the three scenarios
    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_' task '_' num2str(v) '_' storm_name];
            output_elv_name = ['TrElev_5_' task '_' num2str(v) '_' storm_name];
            output_dep_name = ['TrDep_5_' task '_' num2str(v) '_' storm_name];
            output_CD_name = ['TrCD_5_' task '_' num2str(v) '_' 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('A:\Newbury_Modeling\Transect_points')
           % 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')';

                    N = 3093;
                    bv = 0.0027;
                    hv = 0.27;
                    if v == 1 %no vegetation, platform stays the same
                        [WH,El,Dp] = transect_wave_attenuation_novegetation(An,N, bv, ub0, T0, H0, h0, bl0, hv, k0, elv, wave_threshold);
                    elseif v == 2 %no vegetation, platform lowers
                        elv = elv-(abs(elv).*0.1); 
                        [WH,El,Dp] = transect_wave_attenuation_novegetation(An,N, bv, ub0, T0, H0, h0, bl0, hv, k0, elv, wave_threshold);
                    elseif v == 3 %vegtation, platform lowers
                        elv = elv-(abs(elv).*0.1); 
                        [WH,El,Dp,CDa] = transect_wave_attenuation(An,N, bv, ub0, T0, H0, h0, bl0, hv, k0, elv, wave_threshold);
                    elseif v == 4
                        elv = elv+0.25; %(abs(elv).*0.1); 
                        [WH,El,Dp] = transect_wave_attenuation_novegetation(An,N, bv, ub0, T0, H0, h0, bl0, hv, k0, elv, wave_threshold);
                    end
                    R_WA(1:length(WH),j) = WH;
                    R_Elev(1:length(El),j) = El;
                    R_Dep(1:length(Dp),j) = Dp;
                    if v ==3 
                        R_CDa(j) = CDa;
                    end
                    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') ;
            if v ==3 
                save(output_CD_name, 'R_CDa');
            end



        end
    end
end
toc