diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..2a8a2db --- /dev/null +++ b/README.txt @@ -0,0 +1,13 @@ +double click on the file carboCATGUI.m to open Matlab and load this file. +Run the function either by enterting the function name carboCATGUI at the command prompt, +or by clicking run in the editor tab + +The GUI should be reasonably intuitive; click on the Initialise button to set up the model +ready to run, then click on Run CA Model, then on PLot run and you should see the +results from the paramsProcess.txt and the paramsInputValues.txt input files I have set up +in the params folder. You can edit both of these files to run different models. +You can also use the code in the utilities folder to make new initial topography, +intial facies maps, sea-level curves and other inputs you might need to create +useful models + +Any questions, e-mail peter.burgess@liverpool.ac.uk \ No newline at end of file diff --git a/calculateErrorBetweenMarginTrajectory.m b/calculateErrorBetweenMarginTrajectory.m new file mode 100644 index 0000000..8924f1f --- /dev/null +++ b/calculateErrorBetweenMarginTrajectory.m @@ -0,0 +1,207 @@ +clear all +clc + +%%load saved data +strata=load('C:\EquinorProject\Meeting\28Apr\WaveNorth\strata.mat'); +glob.strata=strata.strata; +numberOfLayers=load('C:\EquinorProject\Meeting\28Apr\WaveNorth\numberOfLayers.mat'); +glob.numberOfLayers=numberOfLayers.numberOfLayers; +facies=load('C:\EquinorProject\Meeting\28Apr\WaveNorth\facies.mat'); +glob.facies=facies.facies; +glob.ySize=size(glob.strata,1); +glob.xSize=size(glob.strata,2); +iteration=size(glob.strata,3); +glob.dx=250; %cell size +deltaT=100; %calculate platform margin at every 100 iterations +T=linspace(deltaT,iteration,iteration/deltaT); + +% %%for diagonal cross section where y=x; +% long=min(glob.xSize,glob.ySize); +% gradDiag=zeros(long,iteration); +% PfDiag=zeros(long,iteration); +% Zdiag=zeros(long,iteration); +% for y=1:long-1 +% xPos=y; +% for t=2:iteration +% +% g=glob.strata(y+1,xPos+1,t)-glob.strata(y,xPos,t); +% gradDiag(y,t)=g/25; +% for m=1:glob.numberOfLayers(y,xPos,t) +% if glob.facies{y+1,xPos,t}(1)==10 && glob.facies{y,xPos,t}(m)==1 +% PfDiag(y,t)=1; +% Zdiag(y,t)=glob.strata(y,xPos,t); +% elseif glob.facies{y,xPos,t}(m)==10 && glob.facies{y,xPos,t}(m)==11 && glob.facies{y+1,xPos,t}(1)==1 +% PfDiag(y,t)=1; +% end +% end +% +% +% end +% +% end +% gradDiag(gradDiag>1)=1; +% gradDiag=abs(gradDiag); +% for i=1:length(T) +% t=T(i); +% PDiag(:,i)=gradDiag(:,t).*PfDiag(:,t); +% end +% [Xdiag,Zdiag] = find(PDiag); +% for i=1:size(Xdiag,1) +% Z2diag(i)=glob.strata(Xdiag(i),Xdiag(i),Zdiag(i)*deltaT); +% end +% figure +% subplot(2,1,1) +% scatter(Xdiag,Zdiag,200,'r'); +% ylabel('EMT','FontSize',18,'FontWeight','bold'); +% xlabel('km','FontSize',18,'FontWeight','bold'); +% title(['Diagonal, t= ', num2str(0.001*iteration),'My'],'FontSize',20); +% xticks([1, 0.25*(long+2), 0.5*(long+2), 0.75*(long+2),long+1]); +% xticklabels (round([0 ,long*1.414*glob.dx/4000 ,long*1.414*glob.dx/2000,long*1.414*3*glob.dx/4000,long*1.414*glob.dx/1000])); +% % yticks([0 ,5, 10]); +% % yticklabels ([0 0.5*iteration/1000,iteration/1000]); +% axis([1 long+1 0 iteration/deltaT]) +% ax=gca; +% ax.LineWidth = 0.6; +% ax.FontSize = 24; +% ax.FontWeight = 'bold'; +% +% subplot(2,1,2) +% scatter(Xdiag,Z2diag,200,'r'); +% ylabel('Elevation (m)','FontSize',18,'FontWeight','bold'); +% xlabel('km','FontSize',18,'FontWeight','bold'); +% xticks([1, 0.25*(long+2), 0.5*(long+2), 0.75*(long+2),long+1]); +% xticklabels (round([0 ,long*1.414*glob.dx/4000 ,long*1.414*glob.dx/2000,long*1.414*3*glob.dx/4000,long*1.414*glob.dx/1000])); +% % yticks([0 ,5, 10]); +% % yticklabels ([0 0.5*iteration/1000,iteration/1000]); +% axis([1 long+1 0 1500]) +% ax=gca; +% ax.LineWidth = 0.6; +% ax.FontSize = 24; +% ax.FontWeight = 'bold'; + + +%%for othogonal cross section along y axis +xPos=glob.xSize/2; +gradY=zeros(glob.ySize,iteration); %to store elevation difference between adjacent cells +Pfy=zeros(glob.ySize,iteration); %to store probability of facies transition +Zy=zeros(glob.ySize,iteration); %to store elevation values of platform margins +for y=1:glob.ySize-1 + + for t=2:iteration + + g=glob.strata(y+1,xPos,t)-glob.strata(y,xPos,t); + gradY(y,t)=g/50; + for m=1:glob.numberOfLayers(y,xPos,t) + if glob.facies{y+1,xPos,t}(1)==10 && glob.facies{y,xPos,t}(m)==1 + Pfy(y,t)=1; + Zy(y,t)=glob.strata(y,xPos,t); + elseif glob.facies{y,xPos,t}(m)==10 && glob.facies{y,xPos,t}(m)==11 && glob.facies{y+1,xPos,t}(1)==1 + Pfy(y,t)=1; + end + end + + + end +end +%turn the calculated gradient into a value between 0 and 1 +gradY(gradY>1)=1; +gradY=abs(gradY); +for i=1:length(T) + t=T(i); + Py(:,i)=gradY(:,t).*Pfy(:,t); +end +[X1,Z] = find(Py); +for i=1:size(X1,1) + Z2(i)=glob.strata(X1(i),xPos,Z(i)*deltaT); +end + +%%plot figure in time +figure +subplot('Position',[0.14, 0.55, 0.72 0.35]) +scatter(X1,Z,200,'r'); +ylabel('EMT','FontSize',18,'FontWeight','bold'); +xlabel('y(km)','FontSize',18,'FontWeight','bold'); +title(['t= ', num2str(0.001*iteration),'My'],'FontSize',20); +xticks([1 10 20 30 40 50 60 70 80 90 100]); +xticklabels (0:(glob.dx/100):20); +% yticks([0 ,5, 10]); +% yticklabels ([0 0.5*iteration/1000,iteration/1000]); +axis([1 glob.ySize 0 iteration/deltaT]) +ax=gca; +ax.LineWidth = 0.6; +ax.FontSize = 24; +ax.FontWeight = 'bold'; + +%%for othogonal cross section along x axis +yPos=glob.ySize/2; +gradX=zeros(glob.xSize,iteration); +Pfx=zeros(glob.xSize,iteration); +for x=1:glob.xSize-1 + + for t=2:iteration + g=glob.strata(yPos,x+1,t)-glob.strata(yPos,x,t); + gradX(x,t)=g/50; + for m=1:glob.numberOfLayers(yPos,x,t) + if glob.facies{yPos,x+1,t}(1)==10 && glob.facies{yPos,x,t}(m)==1 + Pfx(x,t)=1; + elseif glob.facies{yPos,x,t}(m)==10 && glob.facies{yPos,x,t}(m)==11 && glob.facies{yPos,x+1,t}(1)==1 + Pfx(x,t)=1; + end + end + + + end +end +gradX(gradX>1)=1; +gradX=abs(gradX); +[row, col] = find(Pfx(:,1:iteration)>0); +for i=1:length(T) + t=T(i); + Px(:,i)=gradX(:,t).*Pfx(:,t); +end + +[X2,Z] = find(Px); +for i=1:size(X2,1) + Z3(i)=glob.strata(yPos,X2(i),Z(i)*deltaT); +end +subplot('Position',[0.14, 0.1, 0.68 0.35]) +scatter(X2,Z,200,'r'); +ylabel('EMT','FontSize',18,'FontWeight','bold'); +xlabel('x(km)','FontSize',18,'FontWeight','bold'); +xticks([1 10 20 30 40 50 60 70 80 90 100]); +xticklabels (0:(glob.dx/100):20); +% yticks([0 ,5, 10]); +% yticklabels ([0 0.5*iteration/1000,iteration/1000]); +axis([1 glob.xSize 0 iteration/deltaT]) +ax=gca; +ax.LineWidth = 0.6; +ax.FontSize = 24; +ax.FontWeight = 'bold'; + + +%%plot cross section in depth +figure +subplot('Position',[0.14, 0.55, 0.72 0.35]) +scatter(X1,Z2,200,'r'); +ylabel('Elevation (m)','FontSize',18,'FontWeight','bold'); +xlabel('y(km)','FontSize',18,'FontWeight','bold'); +title(['t= ', num2str(0.001*iteration),'My'],'FontSize',20); +xticks([1 10 20 30 40 50 60 70 80 90 100]); +xticklabels (0:(glob.dx/100):20); +axis([1 glob.ySize 0 1500]) +ax=gca; +ax.LineWidth = 0.6; +ax.FontSize = 24; +ax.FontWeight = 'bold'; + +subplot('Position',[0.14, 0.1, 0.68 0.35]) +scatter(X2,Z3,200,'r'); +ylabel('Elevation (m)','FontSize',18,'FontWeight','bold'); +xlabel('x(km)','FontSize',18,'FontWeight','bold'); +xticks([1 10 20 30 40 50 60 70 80 90 100]); +xticklabels (0:(glob.dx/100):20); +axis([1 glob.xSize 0 1500]) +ax=gca; +ax.LineWidth = 0.6; +ax.FontSize = 24; +ax.FontWeight = 'bold'; \ No newline at end of file diff --git a/calculateFaciesCAProdDependent.m b/calculateFaciesCAProdDependent.m new file mode 100644 index 0000000..916be8a --- /dev/null +++ b/calculateFaciesCAProdDependent.m @@ -0,0 +1,174 @@ +function [glob] = calculateFaciesCAProdDependent(glob, iteration) +% Calculate the facies cellular automata according to neighbour rules in +% glob.CARules and modify production rate modifier for each cell according to number of neighbours. +% The order in which facies are checked depends on the production value at the depth of the cell + + j = iteration - 1; + k = iteration; + + % Neighbours contains a count of neighbours for each facies at each grid point and is populated by the countAllNeighbours function + neighbours = zeros(glob.ySize, glob.xSize, glob.maxProdFacies); + + [neighbours] = countAllNeighbours(glob, j, neighbours); + glob.faciesProdAdjust = zeros(glob.ySize, glob.xSize); + + tempFaciesProd = zeros(glob.ySize, glob.xSize, glob.maxProdFacies); % need a CA array for each prod facies + + for y = 1 : glob.ySize + for x= 1 : glob.xSize + + oneFacies = glob.facies{y,x,j}(1); + + % Only do anything here if the latest stratal surface is below sea-level, i.e. + % water depth > 0.001 and if the concentration is enough to maintain production + if glob.wd(y,x,k) > 0.001 + + % For a subaerial hiatus, now reflooded because from above wd > 0 + if oneFacies == 9 % 9 is the code for subaerial exposure + checkProcess=strcmp(glob.refloodingRoutine,'pre-exposed'); + if checkProcess==1 + % Get the facies from the pre-exposed + glob.facies{y,x,k}(1) = findPreHiatusFacies(glob, x,y,iteration); % reoccupy with facies from below hiatus + glob.numberOfLayers(y,x,k)=1; + else + % Set facies =0; + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + if glob.facies{y,x,k}(1)> glob.maxProdFacies || glob.facies{y,x,k}(1)==0 + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + + % For cells already containing producing facies + elseif (oneFacies <= glob.maxProdFacies && oneFacies>0) && glob.numberOfLayers(y,x,j)==1 + + % Check if neighbours is less than min for survival CARules(i,2) or + % greater than max for survival CARules(i,3), or if water depth is greater than production cutoff and if so kill + % that facies + thick = glob.carbonateVolMap(y,x) / (glob.ySize*glob.xSize); + checkProcess=strcmp(glob.concentrationRoutine,'on'); + + if (neighbours(y,x,oneFacies) < glob.CARules(oneFacies,2)) ||... + (neighbours(y,x,oneFacies) > glob.CARules(oneFacies,3)) ||... + (glob.wd(y,x,iteration) >= glob.prodRateWDCutOff(oneFacies)) ||... + ( checkProcess == 1 && thick < 0.01) + glob.facies{y,x,k}(1) = 0; + else % The right number of neighbours exists + glob.facies{y,x,k}(1) = oneFacies; + glob.numberOfLayers(y,x,k)=1; + end + + else % Otherwise cell must be empty or contain transported product so see if it can be colonised with a producing facies + % Check again the carbonate concentration + thick = glob.carbonateVolMap(y,x) / (glob.ySize*glob.xSize); + checkProcess=strcmp(glob.concentrationRoutine,'on'); + if checkProcess == 0 || thick > 0.01 + for f = 1:glob.maxProdFacies + % Check if number of neighbours is within range to trigger new + % facies cell, and only allow new cell if the water depth is less the production cut off depth + if (neighbours(y,x,f) >= glob.CARules(f,4)) && (neighbours(y,x,f) <= glob.CARules(f,5)) && (glob.wd(y,x,iteration) < glob.prodRateWDCutOff(f)) + tempFaciesProd(y,x,f) = f; % new facies cell triggered at y,x + end + end + + prodRate = zeros(1,glob.maxProdFacies); + for f = 1:glob.maxProdFacies + if tempFaciesProd(y,x,f) > 0 + prodRate(f) = glob.prodRate(f) * calculateProductionWDAdjustment(glob, x,y, f, iteration); + end + end + + % Find the highest production rate and assign a fraction of 1 at the correct position in the glob.fraction array + [sortedArray,index]=sort(prodRate); + if sum(sortedArray(:))==0 + glob.facies{y,x,k}(1)=0; + else + glob.facies{y,x,k}(1)=index(glob.maxProdFacies); % the last element in the array has the largest production + glob.numberOfLayers(y,x,k)=1; + end + end + end + + % Finally, calculate the production adjustment factor for the current facies distribution + oneFacies = glob.facies{y,x,k}(1); + glob.faciesProdAdjust(y,x)=calculateFaciesProductionAdjustementFactor(glob,y,x,oneFacies,neighbours); + else + glob.facies{y,x,k}(1) = 9; % Set to above sea-level facies because must be at or above sea-level + glob.numberOfLayers(y,x,k)=1; + + end + end + end +end + +function [neighbours] = countAllNeighbours(glob, j, neighbours) +% Function to count the number of cells within radius containing facies 1 +% to maxFacies across the whole grid and store results in neihgbours matrix + + ySize=int16(glob.ySize); + xSize=int16(glob.xSize); + for yco = 1 : ySize + for xco = 1 : xSize + oneFacies = glob.facies{yco,xco,j}(1); + if oneFacies>0 && oneFacies<=glob.maxProdFacies + radius = glob.CARules(oneFacies,1); + else + radius = glob.CARules(1,1); + end + for l = -radius : radius + for m = -radius : radius + + y = yco + l; + x = xco + m; + + checkProcess=strcmp(glob.wrapRoutine,'unwrap'); + if checkProcess==1 + %if near the edge, complete in a mirror-like image the + %neighbours array + if y<1; y=1+(1-y); end + if x<1; x=1+(1-x); end + if y>ySize; y=ySize+(ySize-y); end + if x>xSize; x=xSize+(xSize-x); end + else + %or wrap around the corners + if y<1; y=ySize+1+l; end + if x<1; x=xSize+1+m; end + if y>ySize; y=l; end + if x>xSize; x=m; end + end + %now count the neighbours using the x-y indeces + + %Don't include cell that has over a certain wd + %difference, and penalise the CA with a factor + wdDiff=abs(glob.wd(yco,xco,j)-glob.wd(y,x,j)); + + if wdDiff>glob.BathiLimit + wdDiffFactor=0; + else + wdDiffFactor=(-1/glob.BathiLimit)*wdDiff+1; + end + if wdDiffFactor>0.9; wdDiffFactor=1; end + faciesType=glob.facies{y,x,j}(1); + % Count producing facies as neighbours but do not include the center cell - + % neighbours count should not include itself + if faciesType > 0 && faciesType <= glob.maxProdFacies && not (l == 0 && m == 0) + neighbours(yco,xco,faciesType) = neighbours(yco,xco,faciesType) + (1*wdDiffFactor); + end + + end + end + end + end +end + +function [preHiatusFacies] = findPreHiatusFacies(glob, x,y,iteration) + + k = iteration - 1; + + while k > 0 && glob.facies{y,x,k}(1) == 9 + k = k - 1; + end + + preHiatusFacies = glob.facies{y,x,k}(1); +end diff --git a/calculateFaciesCARotationOrder.m b/calculateFaciesCARotationOrder.m new file mode 100644 index 0000000..83baaeb --- /dev/null +++ b/calculateFaciesCARotationOrder.m @@ -0,0 +1,223 @@ +function [glob] = calculateFaciesCARotationOrder(glob, iteration, order) +% Calculate the facies cellular automata according to neighbour rules in +% glob.CARules and modify production rate modifier for each cell according to number of neighbours +% To avoid bias, facies must be dealt with in different order each time, hence orderArray + + switch glob.maxProdFacies + case 1 + orderArray = [1]; + case 2 + orderArray = [1,2; 2,1]; + case 3 + orderArray = [1,2,3; 2,3,1; 3,1,2]; + otherwise + orderArray = [1,2,3,4; 2,3,1,4; 4,3,2,1;3,4,1,2]; + + % need to add general condition here to deal with 4+ facies cases + % Maybe replace this code with simpler version that just uses an order array but rotates it each iteration?? + end + + j = iteration - 1; + k = iteration; + % Neighbours contains a count of neighbours for each facies at each grid + % point and is populated by the countAllNeighbours function + neighbours = zeros(glob.ySize, glob.xSize, glob.maxProdFacies); + neighbours = countAllNeighbours(glob, j, neighbours); + + glob.faciesProdAdjust = zeros(glob.ySize, glob.xSize); + + for y = 1 : glob.ySize + for x= 1 : glob.xSize + + % only do anything here if the latest stratal surface is below sea-level, i.e. + % water depth > 0.001 + if glob.wd(y,x,iteration) > 0.001 + + % Get the facies currently present at y x for previous iteration + oneFacies = glob.facies{y,x,j}(1); + + if oneFacies == 9 % So a subaerial hiatus, now reflooded because from above wd > 0 + + checkProcess=strcmp(glob.refloodingRoutine,'pre-exposed'); + if checkProcess==1 + %get the facies from the pre-exposed + glob.facies{y,x,k}(1) = findPreHiatusFacies(glob, x,y,iteration); % reoccupy with facies f + glob.numberOfLayers(y,x,k)=1; + else + %set facies =0; + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + if glob.facies{y,x,k}(1)> glob.maxProdFacies || glob.facies{y,x,k}(1)==0 + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + % For cells already containing producing facies (and no + % transported on top) + elseif (oneFacies > 0 && oneFacies <= glob.maxProdFacies) && glob.numberOfLayers(y,x,j)==1 + + % Check if neighbours is less than min for survival CARules(i,2) or + % greater than max for survival CARules(i,3), or if water depth is greater than production cutoff and if so kiil + % that facies or if the concentration of carbonate is too low + thick = glob.carbonateVolMap(y,x) / (glob.ySize*glob.xSize); + checkProcess=strcmp(glob.concentrationRoutine,'on'); + + if (neighbours(y,x,oneFacies) < glob.CARules(oneFacies,2)) ||... + (neighbours(y,x,oneFacies) > glob.CARules(oneFacies,3)) ||... + (glob.wd(y,x,iteration) >= glob.prodRateWDCutOff(oneFacies)) ||... + ( checkProcess == 1 && thick < 0.01) + glob.facies{y,x,k}(1) = 0; + else % The right number of neighbours exists + % Facies persists if wave energy routine is off, or the platform margin is below wave break depth + checkProcess=strcmp(glob.waveRoutine,'on'); + if checkProcess==0 || glob.wd(stats.averagePos(iteration-1,x),x, iteration-1) > (3/0.78) + glob.facies{y,x,k}(1) = oneFacies; + glob.numberOfLayers(y,x,k) = 1; + else% (wave energy is on and the platform margin is above the wave break depth)the facies remains if the right wave energy exists + [modificationFactor] = calculateProductionWaveEnergyAdjustment(glob, x,y, oneFacies); + if modificationFactor == 1 + glob.facies{y,x,k}(1) = oneFacies; + glob.numberOfLayers(y,x,k)=1; + else + glob.facies{y,x,k}(1) = 0; + end + end + end + else % Otherwise cell must be empty or contain transported product so see if it can be colonised with a producing facies + % Note that we do not want iteration count to apply to empty cells, + % hence this else + %check again the carbonate concentration + thick = glob.carbonateVolMap(y,x) / (glob.ySize*glob.xSize); + checkProcess=strcmp(glob.concentrationRoutine,'on'); + if checkProcess == 0 || thick > 0.01 + for m = 1:glob.maxProdFacies % Loop through all the facies + + % Select a facies number from the order array. Remember this makes + % sure that facies occurrence in adjacent cells is checked in a + % different order each time step + p = orderArray(order, m); + + % Check if number of neighbours is within range to trigger new + % facies cell, and only allow new cell if the water depth is less + % the production cut off depth + + if (neighbours(y,x,p) >= glob.CARules(p,4)) && (neighbours(y,x,p) <= glob.CARules(p,5) && glob.wd(y,x,iteration) < glob.prodRateWDCutOff(p)) + %The right number of cells exist + %The facies colonise if wave energy routine is off + checkProcess=strcmp(glob.waveRoutine,'off'); + if checkProcess ==1 + glob.facies{y,x,k}(1) = p; % new facies cell triggered at y,x + glob.numberOfLayers(y,x,k)=1; + else % or if the energy routine is on and the margin is above wave break + if glob.wd(stats.averagePos(iteration-1,x),x,iteration-1) > (3/0.78) + glob.facies{y,x,k}(1) = p; % new facies cell triggered at y,x + glob.numberOfLayers(y,x,k)=1; + else + modificationFactor = calculateProductionWaveEnergyAdjustment(glob, x,y, oneFacies); + if modificationFactor==1 + glob.facies{y,x,k}(1) = oneFacies; + glob.numberOfLayers(y,x,k)=1; + else + glob.facies{y,x,k}(1) = 0; + end + end + end + end + end + end + + order = order + 1; + if order > glob.maxProdFacies + order = 1; + end + end + + % Finally, calculate the production adjustment factor for the current facies distribution + oneFacies = glob.facies{y,x,k}(1); + glob.faciesProdAdjust(y,x)=calculateFaciesProductionAdjustementFactor(glob,y,x,oneFacies,neighbours); + else + glob.facies{y,x,k}(1) = 9; % Set to above sea-level facies because must be at or above sea-level%NB: change from 7 to 9 to support 4 facies + glob.numberOfLayers(y,x,k)=1; + end + end + end +end + +function [neighbours, numberOfNeighbours] = countAllNeighbours(glob, j, neighbours,numberOfNeighbours) +% Function to count the number of cells within radius containing facies 1 +% to maxFacies across the whole grid and store results in neihgbours matrix + + ySize=int16(glob.ySize); + xSize=int16(glob.xSize); + + for yco = 1 : ySize + for xco = 1 : xSize + %get the facies of the centre cell + oneFacies = glob.facies{yco,xco,j}(1); + if oneFacies>0 && oneFacies<=glob.maxProdFacies + %if the centre is occupied by a producing cell, get the radius + %from the file + radius = glob.CARules(oneFacies,1); + else + %if not use two (or one) + radius = glob.CARules(1,1); + end + + for l = -radius : radius + for m = -radius : radius + + y = yco + l; + x = xco + m; + + checkProcess=strcmp(glob.wrapRoutine,'unwrap'); + if checkProcess==1 + %if near the edge, complete in a mirror-like image the + %neighbours array + + if y<1; y=1+(1-y); end + if x<1; x=1+(1-x); end + if y>ySize; y=ySize+(ySize-y); end + if x>xSize; x=xSize+(xSize-x); end + + else + %or wrap around the corners + if y<1; y=ySize+1+l; end + if x<1; x=xSize+1+m; end + if y>ySize; y=l; end + if x>xSize; x=m; end + end + %now count the neighbours using the x-y indeces + + %Don't include cell that has over a certain wd + %difference, and penalise the CA with a factor + wdDiff=abs(glob.wd(yco,xco,j)-glob.wd(y,x,j)); + + if wdDiff>glob.BathiLimit + wdDiffFactor=0; + else + wdDiffFactor=(-1/glob.BathiLimit)*wdDiff+1; + end + if wdDiffFactor>0.9; wdDiffFactor=1; end + faciesType=glob.facies{y,x,j}(1); + % Count producing facies as neighbours but do not include the center cell - + % neighbours count should not include itself + if faciesType > 0 && faciesType <= glob.maxProdFacies && not (l == 0 && m == 0) + neighbours(yco,xco,faciesType) = neighbours(yco,xco,faciesType) + (1*wdDiffFactor); + end + + end + end + end + end +end + +function [preHiatusFacies] = findPreHiatusFacies(glob, x,y,iteration) + + k = iteration - 1; + + while k > 0 && glob.facies{y,x,k}(1) == 9 + k = k - 1; + end + + preHiatusFacies = glob.facies{y,x,k}(1); +end diff --git a/calculateFaciesCAWaveDependent.m b/calculateFaciesCAWaveDependent.m new file mode 100644 index 0000000..68d1618 --- /dev/null +++ b/calculateFaciesCAWaveDependent.m @@ -0,0 +1,205 @@ +function [glob] = calculateFaciesCAWaveDependent(glob, iteration) +% Calculate the facies cellular automata according to neighbour rules in +% glob.CARules and according to the wave energy preferences of each facies + + if strcmp(glob.waveRoutine,'on') == 1 + + % Initialise arrays needed for this function + neighbours = zeros(glob.ySize, glob.xSize, glob.maxProdFacies); + glob.faciesProdAdjust = zeros(glob.ySize, glob.xSize); + possNewFacies = zeros(glob.maxProdFacies, 1); % vector of facies that could colonise an empty cell - see below + + j = iteration - 1; + k = iteration; + + % Neighbours contains a count of neighbours for each facies at each grid point and is populated by the countAllNeighbours function + [neighbours] = countAllNeighbours(glob, j, neighbours); + + for y = 1 : glob.ySize + for x= 1 : glob.xSize + + oneFacies = glob.facies{y,x,j}(1); + + % Only do anything here if the latest stratal surface is below sea-level, i.e. + % water depth > 0.001 and if the concentration is enough to maintain production + if glob.wd(y,x,k) > 0.001 + + % For a subaerial hiatus, now reflooded because from above wd > 0 + if oneFacies == 9 % 9 is the code for subaerial exposure + checkProcess = strcmp(glob.refloodingRoutine,'pre-exposed'); + if checkProcess==1 + %get the facies from the pre-exposed + glob.facies{y,x,k}(1) = findPreHiatusFacies(glob, x,y,iteration); % reoccupy with facies from below hiatus + glob.numberOfLayers(y,x,k)=1; + else + % set facies =0; + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + + if glob.facies{y,x,k}(1)> glob.maxProdFacies || glob.facies{y,x,k}(1) == 0 + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + + % For cells already containing producing facies + elseif (oneFacies <= glob.maxProdFacies && oneFacies>0) && glob.numberOfLayers(y,x,j)==1 + + % Check if neighbours is less than min for survival CARules(i,2) or greater than max for survival CARules(i,3), + % or if water depth is greater than production cutoff, and if so kill that facies + if (neighbours(y,x,oneFacies) < glob.CARules(oneFacies,2)) ||... + (neighbours(y,x,oneFacies) > glob.CARules(oneFacies,3)) ||... + (glob.wd(y,x,iteration) >= glob.prodRateWDCutOff(oneFacies)) + glob.facies{y,x,k}(1) = 0; % kill the facies because wrong neighbour count or too deep + + else % The right number of neighbours exists + + % Facies persists if wave energy is appropriate + if glob.waveEnergy(y,x) >= glob.prodWaveThresholdLow(oneFacies) && ... + glob.waveEnergy(y,x) <= glob.prodWaveThresholdHigh(oneFacies) + glob.facies{y,x,k}(1) = oneFacies; + glob.numberOfLayers(y,x,k) = 1; + else + glob.facies{y,x,k}(1) = 0; + end + end + + else % Otherwise cell must be empty or contain transported product, so see if it can be colonised with a producing facies + possNewFacies = zeros(glob.maxProdFacies, 1); + possNewFaciesCount = 0; + for f = 1:glob.maxProdFacies % loop to create candidate facies cells in the empty grid cell + % Check if number of neighbours is within range to trigger new facies cell, + % and only allow new cell if the water depth is less the production cut off depth + if (neighbours(y,x,f) >= glob.CARules(f,4)) && ... + (neighbours(y,x,f) <= glob.CARules(f,5)) && ... + (glob.wd(y,x,iteration) < glob.prodRateWDCutOff(f)) + possNewFacies(f) = f; % reset from zero for a facies that meets neighbour and water depth rules + possNewFaciesCount = possNewFaciesCount + 1; + end + end + + % More than one candidate facies so calculate which is optimum for wave energy at this x,y + if possNewFaciesCount > 1 + waveOptimumFacies = 0; % default value to set if no optimum wave energy facies found + minDivergence = 100; % higher than maximum possible, since 0<=wave energy<=1 + for f = 1:glob.maxProdFacies + if possNewFacies(f) > 0 + + divergence = abs(glob.prodWaveOptimum(f) - glob.waveEnergy(y,x)); + if divergence < minDivergence + minDivergence = divergence; + waveOptimumFacies = f; + end + end + end + +% if glob.wd(y,x,k) <= glob.prodRateWDCutOff(waveOptimumFacies) % Finally, check the selected facies is good in this water depth + glob.facies{y,x,k}(1) = waveOptimumFacies; + if waveOptimumFacies > 0 + glob.numberOfLayers(y,x,k)=1; % non zero-facies code here means one layer has been produced at y,x,k + end +% end + elseif possNewFaciesCount == 1 + + singlePossNewFacies = max(possNewFacies(:)); + +% if glob.wd(y,x,k) <= glob.prodRateWDCutOff(singlePossNewFacies) && ... + if glob.waveEnergy(y,x) >= glob.prodWaveThresholdLow(singlePossNewFacies) && ... + glob.waveEnergy(y,x) <= glob.prodWaveThresholdHigh(singlePossNewFacies) + + glob.facies{y,x,k}(1) = singlePossNewFacies; + glob.numberOfLayers(y,x,k)=1; + end + end + end + else + glob.facies{y,x,k}(1) = 9; % Set to above sea-level facies because must be at or above sea-level + glob.numberOfLayers(y,x,k)=1; + end + + % Finally, calculate the production adjustment factor for the current facies distribution + oneFacies = glob.facies{y,x,k}(1); + glob.faciesProdAdjust(y,x)=calculateFaciesProductionAdjustementFactor(glob,y,x,oneFacies,neighbours); + end + end + else + fprintf('Need wave energy process selected in process options to use wave energy sensitive CA\n'); + fprintf('WARNING - no CA calculated\n'); + end + + finalFaciesMap = cell2mat(glob.facies(:,:,k)); + f = 1:9; + faciesTotals(f) = sum(finalFaciesMap(:) == f); + fprintf('CA 1:%d 2:%d 3:%d', faciesTotals(1), faciesTotals(2), faciesTotals(3)); +end + +function [neighbours] = countAllNeighbours(glob, j, neighbours) +% Function to count the number of cells within radius containing facies 1 +% to maxFacies across the whole grid and store results in neihgbours matrix + + ySize=int16(glob.ySize); + xSize=int16(glob.xSize); + for yco = 1 : ySize + for xco = 1 : xSize + oneFacies = glob.facies{yco,xco,j}(1); + if oneFacies>0 && oneFacies<=glob.maxProdFacies + radius = glob.CARules(oneFacies,1); + else + radius = glob.CARules(1,1); + end + for l = -radius : radius + for m = -radius : radius + + y = yco + l; + x = xco + m; + + checkProcess=strcmp(glob.wrapRoutine,'unwrap'); + if checkProcess==1 + %if near the edge, complete in a mirror-like image the + %neighbours array + if y<1; y=1+(1-y); end + if x<1; x=1+(1-x); end + if y>ySize; y=ySize+(ySize-y); end + if x>xSize; x=xSize+(xSize-x); end + else + %or wrap around the corners + if y<1; y=ySize+1+l; end + if x<1; x=xSize+1+m; end + if y>ySize; y=l; end + if x>xSize; x=m; end + end + %now count the neighbours using the x-y indeces + + %Don't include cell that has over a certain wd + %difference, and penalise the CA with a factor + wdDiff=abs(glob.wd(yco,xco,j)-glob.wd(y,x,j)); + + if wdDiff>glob.BathiLimit + wdDiffFactor=0; + else + wdDiffFactor=(-1/glob.BathiLimit)*wdDiff+1; + end + if wdDiffFactor>0.9; wdDiffFactor=1; end + faciesType=glob.facies{y,x,j}(1); + % Count producing facies as neighbours but do not include the center cell - + % neighbours count should not include itself + if faciesType > 0 && faciesType <= glob.maxProdFacies && not (l == 0 && m == 0) + neighbours(yco,xco,faciesType) = neighbours(yco,xco,faciesType) + (1*wdDiffFactor); + end + + end + end + end + end +end + +function [preHiatusFacies] = findPreHiatusFacies(glob, x,y,iteration) + + k = iteration - 1; + + while k > 0 && glob.facies{y,x,k}(1) == 9 + k = k - 1; + end + + preHiatusFacies = glob.facies{y,x,k}(1); +end diff --git a/calculateFaciesCAWaveDependentHybrid.m b/calculateFaciesCAWaveDependentHybrid.m new file mode 100644 index 0000000..a30585d --- /dev/null +++ b/calculateFaciesCAWaveDependentHybrid.m @@ -0,0 +1,253 @@ +function [glob] = calculateFaciesCAWaveDependentHybrid(glob, iteration) +% Calculate the facies cellular automata according to neighbour rules in +% glob.CARules and according to the wave energy preferences of each facies + + if strcmp(glob.waveRoutine,'on') == 1 + + % Initialise arrays needed for this function + neighbours = zeros(glob.ySize, glob.xSize, glob.maxProdFacies); + glob.faciesProdAdjust = zeros(glob.ySize, glob.xSize); + j = iteration - 1; + k = iteration; + + glob.faciesProdAdjust = ones(glob.ySize, glob.xSize); + + % Neighbours contains a count of neighbours for each facies at each grid point and is populated by the countAllNeighbours function + [neighbours] = countAllNeighbours(glob, j, neighbours); + + for y = 1 : glob.ySize + for x= 1 : glob.xSize + + % Only do anything here if the latest stratal surface is below sea-level, i.e. + % water depth > 0.001 and if the concentration is enough to maintain production + if glob.wd(y,x,k) > 0.001 + + [glob, newFacies] = coloniseOrContinueWaveDependentFacies(glob, x, y, k); + + if newFacies == 0 + glob = calcWaveInsensitiveCA(glob, x, y, j, k, neighbours); + end + else + glob.facies{y,x,k}(1) = 9; % Set to above sea-level facies because must be at or above sea-level + glob.numberOfLayers(y,x,k)=1; + end + end + end + else + fprintf('Need wave energy process selected in process options to use wave energy sensitive CA\n'); + fprintf('WARNING - no CA calculated\n'); + end + +% finalFaciesMap = cell2mat(glob.facies(:,:,k)); +% f = 1:9; +% faciesTotals(f) = sum(finalFaciesMap(:) == f); +% fprintf('CA 1:%d 2:%d 3:%d ', faciesTotals(1), faciesTotals(2), faciesTotals(3)); +end + +function [neighbours] = countAllNeighbours(glob, j, neighbours) +% Function to count the number of cells within radius containing facies 1 +% to maxFacies across the whole grid and store results in neihgbours matrix + + ySize=int16(glob.ySize); + xSize=int16(glob.xSize); + for yco = 1 : ySize + for xco = 1 : xSize + oneFacies = glob.facies{yco,xco,j}(1); + if oneFacies>0 && oneFacies<=glob.maxProdFacies + radius = glob.CARules(oneFacies,1); + else + radius = glob.CARules(1,1); + end + for l = -radius : radius + for m = -radius : radius + + y = yco + l; + x = xco + m; + + checkProcess=strcmp(glob.wrapRoutine,'unwrap'); + if checkProcess==1 + %if near the edge, complete in a mirror-like image the + %neighbours array + if y<1; y=1+(1-y); end + if x<1; x=1+(1-x); end + if y>ySize; y=ySize+(ySize-y); end + if x>xSize; x=xSize+(xSize-x); end + else + %or wrap around the corners + if y<1; y=ySize+1+l; end + if x<1; x=xSize+1+m; end + if y>ySize; y=l; end + if x>xSize; x=m; end + end + % now count the neighbours using the x-y indices + + % Don't include cell that has over a certain wd + % difference, and penalise the CA with a factor + wdDiff = abs(glob.wd(yco,xco,j)-glob.wd(y,x,j)); + + if wdDiff>glob.BathiLimit + wdDiffFactor=0; + else + wdDiffFactor=(-1/glob.BathiLimit)*wdDiff+1; + end + + if wdDiffFactor>0.9; wdDiffFactor=1; end + faciesType=glob.facies{y,x,j}(1); + + % Count producing facies as neighbours but do not include the center cell - + % neighbours count should not include itself + if faciesType > 0 && faciesType <= glob.maxProdFacies && not (l == 0 && m == 0) + neighbours(yco,xco,faciesType) = neighbours(yco,xco,faciesType) + (1*wdDiffFactor); + end + end + end + end + end +end + +function [preHiatusFacies] = findPreHiatusFacies(glob, x,y,iteration) + + j = iteration - 1; + + while j > 0 && glob.facies{y,x,j}(1) == 9 + j = j - 1; + end + + preHiatusFacies = glob.facies{y,x,j}(1); +end + +function glob = calcWaveInsensitiveCA(glob, x, y, j, k, neighbours) + + oneFacies = glob.facies{y,x,j}(1); % Get the previously deposited facies at x,y + + % For a subaerial hiatus, now reflooded because from above wd > 0 + if oneFacies == 9 % 9 is the code for subaerial exposure + checkProcess = strcmp(glob.refloodingRoutine,'pre-exposed'); + if checkProcess==1 + %get the facies from the pre-exposed + glob.facies{y,x,k}(1) = findPreHiatusFacies(glob, x,y,k); % reoccupy with facies from below hiatus + glob.numberOfLayers(y,x,k)=1; + else + % set facies =0; + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + + if glob.facies{y,x,k}(1)> glob.maxProdFacies || glob.facies{y,x,k}(1) == 0 + glob.facies{y,x,k}(1) = 0; + glob.numberOfLayers(y,x,k)=0; + end + + % For cells already containing producing facies + elseif (oneFacies <= glob.maxProdFacies && oneFacies > 0) && ... + glob.prodWaveThresholdHigh(oneFacies) - glob.prodWaveThresholdLow(oneFacies) == 1.0 && ... + glob.numberOfLayers(y,x,j)==1 + + % Check if neighbours is less than min for survival CARules(i,2) or greater than max for survival CARules(i,3), + % or if water depth is greater than production cutoff, and if so kill that facies + if (neighbours(y,x,oneFacies) < glob.CARules(oneFacies,2)) ||... + (neighbours(y,x,oneFacies) > glob.CARules(oneFacies,3)) ||... + (glob.wd(y,x,k) >= glob.prodRateWDCutOff(oneFacies)) + glob.facies{y,x,k}(1) = 0; % kill the facies because wrong neighbour count or too deep + + else % The right number of neighbours exists so facies persists + + glob.facies{y,x,k}(1) = oneFacies; + glob.numberOfLayers(y,x,k) = 1; + end + else % Otherwise cell must be empty or contain transported product, so see if it can be colonised with a producing facies + newFaciesProd = zeros(1,glob.maxProdFacies); + for f = 1:glob.maxProdFacies + % Check if number of neighbours is within range to trigger new + % facies cell, and only allow potential new faices to colonise the cell if the water depth is less the production cut off depth + if (neighbours(y,x,f) >= glob.CARules(f,4)) && ... + (neighbours(y,x,f) <= glob.CARules(f,5)) && ... + (glob.wd(y,x,k) < glob.prodRateWDCutOff(f)) && ... + glob.waveEnergy(y,x) >= glob.prodWaveThresholdLow(f) && ... + glob.waveEnergy(y,x) <= glob.prodWaveThresholdHigh(f) + + newFaciesProd(f) = f; % new facies cell triggered at y,x + end + end + + % calculate the water-depth adjusted production rate for each facies and store in prodRate vector + prodRate = zeros(1,glob.maxProdFacies); + for f = 1:glob.maxProdFacies + if newFaciesProd(f) > 0 + prodRate(f) = glob.prodRate(f) * calculateProductionWDAdjustment(glob, x,y, f, k); + end + end + + maxProd = max(prodRate); % get the maximum prod rate and it's index, the latter being the max prod facies code + + if maxProd == 0 % No facies producing at this x,y water depth so cell is empty + glob.facies{y,x,k}(1) = 0; + else + maxProdFacies = find(prodRate == maxProd); + + if length(maxProdFacies) == 1 % only one facies at the maximum prod rate, so assign this facies to cell x,y + glob.facies{y,x,k}(1) = max(prodRate); + else + maxProdFacies = find(prodRate == maxProd); + chooseFacies = randi(length(maxProdFacies),1); + glob.facies{y,x,k}(1) = maxProdFacies(chooseFacies); +% fprintf('\n%d %d new facies %d\n', x,y, maxProdFacies(chooseFacies)); + end + end +% % Find the highest production rate and assign a fraction of 1 at the correct position in the glob.fraction array +% [sortedProdArray,index] = sort(prodRate); +% if sum(sortedProdArray(:)) == 0 +% glob.facies{y,x,k}(1) = 0; +% else +% glob.facies{y,x,k}(1) = index(glob.maxProdFacies); % the last element in the array has the largest production +% glob.numberOfLayers(y,x,k)=1; +% end + end +end + +function [glob, newFacies] = coloniseOrContinueWaveDependentFacies(glob, x, y, k) + + possNewFacies = zeros(glob.maxProdFacies, 1); + possNewFaciesCount = 0; + for f = 1:glob.maxProdFacies % loop through all the producing facies that could be wave-energy dependent + + % First check if facies f is wave-energy sensitive + if glob.prodWaveThresholdHigh(f)- glob.prodWaveThresholdLow(f) < 0.9999 % 0.9999 to allow for rounding error + + % it is wave sensitive, so check that this xy position is + % within the required wave energy and water depth ranges + if glob.waveEnergy(y,x) >= glob.prodWaveThresholdLow(f) && ... % So check if xy pos meets it's WD and wave energy criteria + glob.waveEnergy(y,x) <= glob.prodWaveThresholdHigh(f) && ... + glob.wd(y,x,k) <= glob.prodRateWDCutOff(f) + + possNewFacies(f) = f; % reset from zero for a facies that meets neighbour and water depth rules + possNewFaciesCount = possNewFaciesCount + 1; + end + end + end + + % More than one candidate facies, so calculate which is optimum for the wave energy at this x,y + if possNewFaciesCount > 0 + + % More than one wave sensitive facies could go in this xy cell, so + % select which based on which has an optimum (middle of the range) + % wave energy sensitivity closest to the wave energy at xy + waveOptimumFacies = 0; % default value to set if no optimum wave energy facies found + minDivergence = 100; % higher than maximum possible, since 0<=wave energy<=1 + for f = 1:glob.maxProdFacies + if possNewFacies(f) > 0 + divergence = abs(glob.prodWaveOptimum(f) - glob.waveEnergy(y,x)); + if divergence < minDivergence + minDivergence = divergence; + waveOptimumFacies = f; + end + end + end + + newFacies = waveOptimumFacies; + glob.facies{y,x,k}(1) = waveOptimumFacies; + glob.numberOfLayers(y,x,k) = 1; % non zero-facies code here means one layer has been produced at y,x,k + else + newFacies = 0; + end +end \ No newline at end of file diff --git a/calculateFaciesProductionAdjustementFactor.m b/calculateFaciesProductionAdjustementFactor.m new file mode 100644 index 0000000..f7dc727 --- /dev/null +++ b/calculateFaciesProductionAdjustementFactor.m @@ -0,0 +1,32 @@ +function [faciesProdAdjust]=calculateFaciesProductionAdjustementFactor(glob,y,x,oneFacies,neighbours) + +% Production adjustment is calculated based on +% the number of neighbours for this cell, which has already been +% calculated at the start of this function. +% Note calculation of prod adjust needs to be done regardless of +% iteration count at y x because surrounding cells may have +% iterated and their count changed + +if oneFacies > 0 && oneFacies <= glob.maxProdFacies + + oneNeighbours = neighbours(y,x,oneFacies); + + + if oneNeighbours < glob.prodScaleMin(oneFacies) % So fewer than minimum neighbours + faciesProdAdjust = 0.0; + + elseif oneNeighbours <= glob.prodScaleOptimum(oneFacies) % More than min, fewer than optimum + + faciesProdAdjust = (oneNeighbours-(glob.prodScaleMin(oneFacies)-1))/(glob.prodScaleOptimum(oneFacies)-(glob.prodScaleMin(oneFacies)-1)); + + elseif oneNeighbours <= glob.prodScaleMax(oneFacies) % More than optimum, fewer than max + + faciesProdAdjust = ((glob.prodScaleMax(oneFacies)+1)-oneNeighbours)/((glob.prodScaleMax(oneFacies)+1)-glob.prodScaleOptimum(oneFacies)); + + else % More than maximum number of neighbours + faciesProdAdjust = 0.0; + end +else + faciesProdAdjust=0.0; +end +end \ No newline at end of file diff --git a/calculateFaciesWaveDependentSimple.m b/calculateFaciesWaveDependentSimple.m new file mode 100644 index 0000000..57e31a2 --- /dev/null +++ b/calculateFaciesWaveDependentSimple.m @@ -0,0 +1,29 @@ +function [glob] = calculateFaciesWaveDependentSimple(glob, iteration) +% Calculate the facies cellular automata according to neighbour rules in +% glob.CARules and according to the wave energy preferences of each facies + + if strcmp(glob.waveRoutine,'on') == 1 + + glob.faciesProdAdjust = ones(glob.ySize, glob.xSize); + + for y = 1 : glob.ySize - 1 % Because of y+1 in glob.waveEnergy reference below + for x= 1 : glob.xSize + + if glob.wd(y,x,iteration) > 0.001 + f = 1:glob.maxProdFacies; + waveEnergyImbalance = abs(glob.prodWaveOptimum(f) - glob.waveEnergy(y+1,x)); % y+1 because waves from cell y+1 break in cell y + [~, newFacies] = min(waveEnergyImbalance); + + if (glob.wd(y,x,iteration) < glob.prodRateWDCutOff(newFacies)) + glob.facies{y,x,iteration}(1) = newFacies; + glob.numberOfLayers(y,x,iteration)=1; + else + glob.facies{y,x,iteration}(1) = 0; + end + else + glob.facies{y,x,iteration}(1) = 0; + end + end + end + end +end diff --git a/calculatePlatformMarginTrajectory.m b/calculatePlatformMarginTrajectory.m new file mode 100644 index 0000000..56a0d28 --- /dev/null +++ b/calculatePlatformMarginTrajectory.m @@ -0,0 +1,42 @@ +function platMarginTraj = calculatePlatformMarginTrajectory(glob, iteration, xPos) + + platMarginTraj = zeros(1, iteration); + rateOfChangeOfGradZ = zeros(glob.ySize, iteration); + marginFaciesTransition = zeros(glob.ySize, iteration); + pMap = zeros(glob.ySize, iteration); + waveSensitiveFacies = 1; + + for t=2:iteration + for y = 2:glob.ySize-1 + + % Calculate the rate of change of topographic gradient along the depositional surface + zDiffLeft = glob.strata(y-1,xPos,t) - glob.strata(y,xPos,t); + zDiffRight = glob.strata(y,xPos,t) - glob.strata(y+1,xPos,t); + rateOfChangeOfGradZ(y,t) = (zDiffLeft - zDiffRight) / glob.dx; + + faciesCount = zeros(1,glob.maxFacies); + m = 1:glob.numberOfLayers(y,xPos,t); + oneFacies = glob.facies{y,xPos,t}(m); + faciesCount(oneFacies) = faciesCount(oneFacies) + 1; + [~,mostFreqFacies] = max(faciesCount); + if mostFreqFacies == waveSensitiveFacies + marginFaciesTransition(y,t) = 1; + end + end + end + + rateOfChangeOfGradZ = rateOfChangeOfGradZ / max(max(rateOfChangeOfGradZ)); + pMap = abs(rateOfChangeOfGradZ) .* marginFaciesTransition; + [~,platMarginTraj] = max(pMap); +% platMargintraj = platMarginTraj ; % Scale to correct y coordinate scale + + figure; + yGridVector = (1:glob.ySize) .* glob.dx; + zGridVector = (1:iteration) .* glob.deltaT; + handle = pcolor(zGridVector, yGridVector, pMap); + set(handle, 'EdgeColor', 'none'); + xlabel('Elapsed model time (iteration) (My)'); + ylabel('Y distance (m)'); + hold on + scatter((1:iteration).*glob.deltaT, platMarginTraj.* glob.dx, 'MarkerEdgeColor',[0 .5 .5], 'MarkerFaceColor',[1 .3 0]); +end diff --git a/calculateProduction.m b/calculateProduction.m new file mode 100644 index 0000000..fd6e938 --- /dev/null +++ b/calculateProduction.m @@ -0,0 +1,104 @@ +function glob = calculateProduction(glob, iteration) + + onefacies = uint8(0); + prod = 0.0; + glob.initialTransportThickMap = zeros(glob.ySize, glob.xSize); + availableForTransport = 0.0; + totalProd=0.0; + faciesCellCount = zeros(1,glob.maxProdFacies); + + % Calculate carbonate production for each point on the grid. + % NB could be optimized of prod was put into an xy map then used outside the loop to update all these other xy maps + for y=1:glob.ySize + for x=1:glob.xSize + + oneFacies = glob.facies{y,x,iteration}(1); + + + % Only need to calculate production for occupied cells i.e. 0 < facies < 9 (9 is subaerial exposure) + if oneFacies > 0 && oneFacies <= glob.maxProdFacies && glob.wd(y,x,iteration) > 0 + + faciesCellCount(oneFacies) = faciesCellCount(oneFacies) + 1; % Add the facies cell to the sum count for that facies + + % Bosscher and Schlager production variation with water depth + if (glob.wd(y,x,iteration) > 0.0) + + if glob.surfaceLight(oneFacies) > 500 + % use this if you want to use Schlager's production profile + glob.prodDepthAdjust(y,x) = tanh((glob.surfaceLight(oneFacies) * exp(-glob.extinctionCoeff(oneFacies) * glob.wd(y,x,iteration)))/ glob.saturatingLight(oneFacies)); + else + % use this if you want to use new production profile + glob.prodDepthAdjust(y,x) = (1./ (1+((glob.wd(y,x,iteration)-glob.profCentre(oneFacies))./glob.profWidth(oneFacies)).^(2.*glob.profSlope(oneFacies))) ); + end + else + glob.prodDepthAdjust(y,x) = 0; + end + + % Set production thickness to depth and neighbour-adjusted + % thickness for this facies, taking into account the + % varoius factors that control production + % Note that faciesProdAdjust is calculated in calculateFaciesCA + prod = glob.prodRate(oneFacies) * glob.prodDepthAdjust(y,x) * glob.faciesProdAdjust(y,x) * glob.carbProdCurve(iteration); + prod = prod * glob.productionBySiliciclasticMap(y,x,oneFacies); + + if prod > glob.wd(y,x,iteration) % if production > accommodation set prod=accommodation to avoid build above SL + prod = glob.wd(y,x,iteration); + end + + %check against the dissolved carbonate volume + checkProcess=strcmp(glob.concentrationRoutine,'on'); + if checkProcess==1 + volp=prod*(glob.dx^2); %in metres + if volp>glob.carbonateVolMap(y,x) + volp = glob.carbonateVolMap(y,x); + glob.carbonateVolMap(y,x) = 0; + else + glob.carbonateVolMap(y,x) = glob.carbonateVolMap(y,x)-volp; + end + prod = volp/(glob.dx^2); + end + + % Decrease WD by amount of production + glob.wd(y,x,iteration) = glob.wd(y,x,iteration) - prod; + + % Calculate the transportable thickness at this point as a proportion of production. + checkProcess=strcmp(glob.waveRoutine,'on'); + if checkProcess ==0 + availableForTransport = prod * glob.transportFraction(oneFacies); + else + availableForTransport = prod * calculateTransportFraction(y,x,oneFacies,glob,iteration); + end + glob.initialTransportThickMap(y,x) = availableForTransport; + + % Record the production as thickness in the strata array + glob.strata(y,x,iteration) = glob.strata(y,x,iteration-1) + prod; + glob.thickness{y,x,iteration}(1) = prod; + + + else % No deposition so record zero thickness + if glob.facies{y,x,iteration}(1)~=8 + glob.prodDepthAdjust(y,x) = 0; + prod = 0; + glob.strata(y,x,iteration) = glob.strata(y,x,iteration-1); + glob.thickness{y,x,iteration}(1) = 0.0; + end + end + +% if x==round(glob.xSize / 2) && y == round(glob.ySize / 2) +% fprintf('SL %3.2f WD:%3.2f @%d,%d Facies %d Prod %3.2f ', glob.SL(iteration), glob.wd(y,x,iteration), x,y, oneFacies, prod); +% end + + totalProd = prod + totalProd; + end + end + + fprintf('SL %3.2f WD:%3.2f ', glob.SL(iteration), glob.wd(y,x,iteration)); + for f = 1:glob.maxProdFacies + fprintf('%d:%d ', f, faciesCellCount(f)); + end + fprintf('Tot accum. %3.2f ', totalProd); +end + + + + diff --git a/calculateProductionWDAdjustment.m b/calculateProductionWDAdjustment.m new file mode 100644 index 0000000..491b745 --- /dev/null +++ b/calculateProductionWDAdjustment.m @@ -0,0 +1,8 @@ +function [prodAdjust] = calculateProductionWDAdjustment(glob, x,y, oneFacies, iteration) + + if glob.surfaceLight(oneFacies)>500 + prodAdjust = tanh((glob.surfaceLight(oneFacies) * exp(-glob.extinctionCoeff(oneFacies) * glob.wd(y,x,iteration)))/ glob.saturatingLight(oneFacies)); + else + prodAdjust = (1./ (1+((glob.wd(y,x,iteration)-glob.profCentre(oneFacies))./glob.profWidth(oneFacies)).^(2.*glob.profSlope(oneFacies))) ); + end +end \ No newline at end of file diff --git a/calculateProductionWaveEnergyAdjustment.m b/calculateProductionWaveEnergyAdjustment.m new file mode 100644 index 0000000..db65261 --- /dev/null +++ b/calculateProductionWaveEnergyAdjustment.m @@ -0,0 +1,9 @@ +function modificationFactor = calculateProductionWaveEnergyAdjustment(glob, x,y, faciesCode) + + modificationFactor = 0; + if faciesCode > 0 && faciesCode <= glob.maxProdFacies + if glob.wave(y,x) >= glob.prodWaveThresholdLow(faciesCode) && glob.wave(y,x) <= glob.prodWaveThresholdHigh(faciesCode) + modificationFactor = 1; + end + end +end \ No newline at end of file diff --git a/calculateSoilDeposition.m b/calculateSoilDeposition.m new file mode 100644 index 0000000..0ac390c --- /dev/null +++ b/calculateSoilDeposition.m @@ -0,0 +1,78 @@ +function glob = calculateSoilDeposition (glob,iteration) +% deposit soil and diffuse it + +rate=0; %m My-1 +thicknessToDeposit=rate.*glob.deltaT; +soilThickness=zeros(glob.ySize,glob.xSize); +wdMap=glob.SL(iteration)-glob.strata(:,:,iteration); +%calculate +soilThickness(wdMap<=0)=thicknessToDeposit; +newThickness=zeros(glob.ySize,glob.xSize); +%add to topo + + + +%diffuse +k=500; %in m^2 yr-1 +%scale + +loopT=glob.deltaT*1000000; + +testDiffusion=((glob.dx)^2)/(3*loopT); %in m^2 yr-1 + +%check if the diff coeff is small enough +while k>=testDiffusion +loopT=loopT/10; %decrease deltaT +testDiffusion=((glob.dx)^2)/(3*loopT); +end + +totalLoops=(glob.deltaT*1000000)/loopT; +k=k*loopT; +%get the elevation after production and deposition od suspended sediment +elev = glob.strata(:,:,iteration)+soilThickness; + +originalBasement = glob.strata(:,:,iteration); + +basement=originalBasement; + +for l=1:totalLoops +%Calculate the second derivative in each dimension +dzdx2 = circshift(elev,[0,-1]) + circshift(elev,[0,1]) - (2.*elev); +dzdy2 = circshift(elev,[-1,0]) + circshift(elev,[1,0]) - (2.*elev); + +%Diffuse + +deltaHx = (k.*dzdx2) ./ ((2*glob.dx)^2); +deltaHy = (k.*dzdy2) ./ ((2*glob.dx)^2); + +%create new elevation by adding each dimension +elevAfterDiffusion = (elev+deltaHy+deltaHx); +elev=elevAfterDiffusion; +end + + +thickMap=elevAfterDiffusion-originalBasement; +thickMap(thickMap<0)=0; +%correct for exposed cells + +for y=1:glob.ySize + for x=1:glob.xSize + + if glob.facies{y,x,iteration}(1)==9; + if thickMap(y,x)>0 + if glob.numberOfLayers(y,x,iteration)>1 + glob.thickness{y,x,iteration}(glob.numberOfLayers(y,x,iteration)+1)=thickMap(y,x); + glob.strata(y,x,iteration)=glob.strata(y,x,iteration)+thickMap(y,x); + glob.facies{y,x,iteration}(glob.numberOfLayers(y,x,iteration)+1)=9; + else + glob.thickness{y,x,iteration}(1)=thickMap(y,x); + glob.strata(y,x,iteration)=glob.strata(y,x,iteration)+thickMap(y,x); + glob.facies{y,x,iteration}(1)=9; + end + end + + end + end +end + +end diff --git a/calculateSubsidence.m b/calculateSubsidence.m new file mode 100644 index 0000000..9a492a5 --- /dev/null +++ b/calculateSubsidence.m @@ -0,0 +1,12 @@ +function [glob] = calculateSubsidence(glob, iteration) +% Calculate subsidence over the whole model grid for all strata up to current iteration +% Assumes that subs.subRateMap is the correct maxX and maxY size and contains the subsidence rate appropriate for the model timestep + + % Implicit loop to calculate the subsidence rate to all previous layers of strata, up to current iteration + k = 1:iteration; + glob.strata(:,:,k) = glob.strata(:,:,k) - glob.subRateMap; +end + + + + diff --git a/calculateTransport.m b/calculateTransport.m new file mode 100644 index 0000000..8a51d2e --- /dev/null +++ b/calculateTransport.m @@ -0,0 +1,284 @@ +function [glob] = calculateTransport(glob, iteration) + + topog = zeros(glob.ySize, glob.xSize); + flowThick = 0; + prodFacies = uint8(0); + oneTransThickMap = num2cell(zeros(glob.ySize, glob.xSize)); % Map of thickness of deposited transported sediment for current iteration + oneTransFaciesMap = num2cell(uint16(zeros(glob.ySize, glob.xSize))); % Map of facies deposited by transport for current iteration + transFacies = uint16(0); + flowCount = uint16(0); + flowRoute = zeros(glob.ySize, glob.xSize); + destCells = uint16(0); + noRoute = uint16(0); + foundDeeperCell = false; + dummy1 = uint16(0); + dummy2 = uint16(0); + + topog = glob.strata(:,:,iteration); + + % Calculate transport from each point on the grid + yArray=1:glob.ySize; + xArray=1:glob.xSize; + yArrayRand = yArray(randperm(length(yArray))); + xArrayRand = xArray(randperm(length(xArray))); + for i=1:glob.ySize + for j=1:glob.xSize + + y=yArrayRand(i); + x=xArrayRand(j); + + flowThick = glob.initialTransportThickMap(y,x); % Trans thick map contains the transportable + %thicknesses calculated in the production routines, derived from the glob.erosionPerc value + prodFacies = glob.facies{y,x,iteration}(1); % Get the in-situ produced facies for this cell + + % Check if there is a deeper cell adjacent to source cell before triggering a + % flow. This is necessary because triggering a flow causes sediment to be + % removed from the source cell - don't want this to happen if flow is not + % going to happen + [foundDeeperCell, deepestX, deepestY] = findDeepestNeighbourCell(glob,topog,x,y); + flowGradient = topog(y,x) - topog(deepestY, deepestX); + + % Only call the flow calculator if... + % pTransportOverride < glob.transportContinueProbability OR all the following are true: + % Source cell contrains primary producer facies 1,2, or 3 that are below sea-level + % There is a deeper adjacent cell + % gradient to adjacent deeper cell and flow thickness both exceed set thresholds + if (foundDeeperCell == true ... + && prodFacies <= glob.maxProdFacies && prodFacies > 0 && glob.wd(y,x,iteration) >= 0 ... + && flowGradient > glob.transportGradient(prodFacies) && flowThick > 0.001) + + transFacies = glob.transportProductFacies(prodFacies); % Get the transport facies produced from the in-situ facies + + % Remove the transportable thickness from the relevant array records + glob.strata(y,x,iteration) = glob.strata(y,x,iteration) - flowThick; + glob.thickness{y,x,iteration}(1) = glob.thickness{y,x,iteration}(1) - flowThick; + + [topog, glob] = calcOneTransportEvent(glob, iteration, topog, y, x, flowThick, prodFacies, transFacies, noRoute); + flowCount = flowCount + 1; + end + end + end + + totalTransported = 0.0; + + for y=1:glob.ySize + for x=1:glob.xSize + + checkProcess=strcmp(glob.siliciclasticsRoutine,'on'); + if checkProcess==1 + %calculate the amount of siliciclastic + kl=find(glob.facies{y,x,iteration}(:)==8); + if isempty(kl) + clast = 0; + else + clast = sum(glob.thickness{y,x,iteration}(kl)); + end + else + clast =0; + end + % Increase deposited thickness by amount of deposition in all the transported facies at y,x + % deposited thickness is the amount of deposited without the siliciclastic + if glob.facies{y,x,iteration}(1) <= glob.maxProdFacies + oneThickness=sum(glob.thickness{y,x,iteration}(:))-sum(glob.thickness{y,x,iteration}(1))-clast; + else + oneThickness=sum(glob.thickness{y,x,iteration}(:))-clast; + end + + % Increase deposited thickness by amount of deposition in all the transported facies at y,x + glob.strata(y,x,iteration) = glob.strata(y,x,iteration) + oneThickness; + totalTransported = totalTransported + oneThickness; + + % Decrease WD by amount of deposition + glob.wd(y,x,iteration) = glob.wd(y,x,iteration) - oneThickness; + + if oneThickness>0; destCells=destCells+1; end + end + end + + % Record transported sediment as thickness in the strata array and as facies in the + % faciesTrans array + %glob.faciesTrans(:,:,iteration) = oneTransFaciesMap; + %glob.faciesTransThick(:,:,iteration) = oneTransThickMap; + + fprintf('Trans %3.2f in %d flows to %d cells. Failed %d', totalTransported, flowCount, destCells, noRoute); +end + +function [topog, glob] = calcOneTransportEvent(glob, iteration, topog, startY, startX, flowThick, prodFacies, transFacies, noRoute) +% Calculates the transport from cell startX startY and records the deposited thickness in oneTransThickMap + + % Initialize variables + flowDone = false; + flowRoute = zeros(glob.ySize, glob.xSize); + flowLength = uint16(0); % length of flow in number of grid cells + numFaciesInCell = uint8(0); + flowRoute(startY, startX) = flowLength; + flowGradient = 0; + + y = startY; + x = startX; + + while flowDone == false + + [foundDeeperCell, deepestX, deepestY] = findDeepestNeighbourCell(glob,topog, x, y); + flowGradient = (topog(y,x) - topog(deepestY, deepestX)) / glob.dx; % dx is in m too, so gradient is dimensionaless ratio + + % pTransportOverride is a random number which if less than glob.transportContinueProbability for each facies + % forces transport to continue, even if there is no deeper neighbour cell + pTransportOverride = rand(1); + + % Lower cells found adjacent to the current flow cell, and gradient > minimum threshold, so deposit and prepare to carry on flow + + % TRANSPORT to deeper cell and beyond, with possibility of some deposition also + if foundDeeperCell == true || pTransportOverride < glob.transContinueProb(prodFacies) + if flowGradient > glob.transportGradient(prodFacies) || pTransportOverride < glob.transContinueProb(prodFacies) % So still above deposition threshold + flowLength = flowLength + 1; + flowRoute(deepestY, deepestX) = flowLength; + x = deepestX; % update flow coorindates to this adjacent deeper cell + y = deepestY; + flowDone = false; % Continue flow, because there may be yet deeper cells adjacent to this one... + else + % Calculate the thickness to deposit based on a proportional decay of + % thickness and a minimum cutoff + if flowThick < 0.01 % gradient too small and flow too thin + depositHere = flowThick; + flowThick = 0; + flowDone = true; + else %gradient too small and thick flow + depositHere = flowThick * 0.5; + flowThick = flowThick - depositHere; + end + + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = transFacies; + glob.thickness{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = depositHere; % Nb this could lead to deposition to above sealevel + topog(y, x) = topog(y, x) + depositHere; + end + + % DEPOSIT + % Else no lower cell found, so deposit all flow at the last deepest cell x y unless it is the start cell + else if not(x == startX && y == startY) + availableSpace = glob.SL(iteration)-topog(y,x); + if flowThick <= availableSpace + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = transFacies; + glob.thickness{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = flowThick; + topog(y, x) = topog(y, x) + flowThick; % NB this means flows will interact with deposits of previous flow + + flowDone = true; % End the flow because no deeper cell found or below gradient threshold so deposit in this low + else + %NB add this on 26/5 to avoid building over sea level - estani & george + depositHere = availableSpace; + flowThick = flowThick - depositHere; + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = transFacies; + glob.thickness{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = depositHere; + topog(y, x) = topog(y, x) + depositHere; + + xInc = [0 -1 1 -1 1 -1 1 0 ]; + yInc = [-1 1 -1 -1 0 0 1 1 ]; + a=1; + checkSurr=0; + + while flowThick>0 + ny=y+yInc(a) ; + if ny<1 ; ny=1; end + if ny>glob.ySize; ny=glob.ySize; end + nx=x+xInc(a); + if nx<1 ; nx=1; end + if nx>glob.xSize; nx=glob.xSize; end + availableSpace=glob.SL(iteration)-topog(ny,nx); + if availableSpace>flowThick + depositHere=flowThick; + else + depositHere=availableSpace; + end + flowThick=flowThick-depositHere; + + if depositHere > 0 + glob.numberOfLayers(ny,nx,iteration)=glob.numberOfLayers(ny,nx,iteration)+1; + glob.facies{ny,nx,iteration}(glob.numberOfLayers(ny,nx,iteration)) = transFacies; + glob.thickness{ny,nx,iteration}(glob.numberOfLayers(ny,nx,iteration)) = depositHere; + topog(ny,nx) = topog(ny,nx) + depositHere; + end + + if availableSpace<=0 + checkSurr=1; + end + + a=a+1; + if a>8 + if checkSurr==1 + flowThick = 0; + else + y=y+1; + a=1; + if y==glob.ySize + flowThick=0; + end + end + end + + end + flowDone=true; + end + else + flowDone = true; % Should only happen if deepest cell is start cell + noRoute = noRoute + 1; + end + end + end +end + +function [foundDeeperCell, deepestX, deepestY] = findDeepestNeighbourCell(glob,topog,x,y) + + foundDeeperCell = false; + deepestCellHeight = topog(y,x); + deepestY = y; % These will be the returned values if no deeper cell is found. Important because used in gradient calc above + deepestX = x; + + xInc = [0 0 1 -1 -1 1 -1 1 ]; + yInc = [-1 1 0 0 1 -1 -1 1 ]; + dummy = 0;% check if out of boundaries is the deepest cell + + for k = 1:8 % Loop through the adjacent cells + + % Calculate xco and yco from the value of x y modifed by the increments in + % xInc and yInc. Should give a coord for each of the adjacent 8 cells + xWrap = x + xInc(k); + yWrap = y + yInc(k); + + if strcmp(glob.wrapRoutine,'unwrap') == 1 + % if on the grid edge, inplement closed boundary boundary + if xWrap < 1; xWrap=1; end; + if xWrap > glob.xSize; xWrap=glob.xSize; end; + if yWrap < 1; yWrap = 1; end; + if yWrap > glob.ySize; yWrap = glob.ySize; end; + else + % if on the grid edge, implement wrapping-mirrored-reflecting etc boundary + if xWrap < 1; xWrap=glob.xSize; end; + if xWrap > glob.xSize; xWrap=1; end; + if yWrap < 1; yWrap = glob.ySize; end; + if yWrap > glob.ySize; yWrap = 1; end; + end + + % Check the current grid cell in the neighbour loop and see if it is the deepest found so far and not already visited + % Cell is deeper + if topog(yWrap, xWrap) < deepestCellHeight + foundDeeperCell = true; + deepestCellHeight = topog(yWrap, xWrap); % added 5 Jan 2014 - elimnates depositional ridges bug on the slope + deepestX = xWrap; + deepestY = yWrap; + + % two (or more(adjacent) have the same depth. Move to the closest, accounting for diagonals + elseif topog(yWrap, xWrap) == deepestCellHeight + foundDeeperCell = true; + distDeep = sqrt(double((y-deepestY)^2 + (x-deepestX)^2)); + distNew = sqrt(double((y-yWrap)^2 + (x-xWrap)^2)); + if distNew <= distDeep + deepestCellHeight = topog(yWrap, xWrap); % added 5 Jan 2014 - elimnates depositional ridges bug on the slope + deepestX = xWrap; + deepestY = yWrap; + end + end + end +end diff --git a/calculateTransportCrossPlatform.m b/calculateTransportCrossPlatform.m new file mode 100644 index 0000000..afc8df8 --- /dev/null +++ b/calculateTransportCrossPlatform.m @@ -0,0 +1,539 @@ +function [glob] = calculateTransportCrossPlatform(glob, iteration) + +topog = zeros(glob.ySize, glob.xSize); +flowThick = 0; +prodFacies = uint8(0); +oneTransThickMap = num2cell(zeros(glob.ySize, glob.xSize)); % Map of thickness of deposited transported sediment for current iteration +oneTransFaciesMap = num2cell(uint16(zeros(glob.ySize, glob.xSize))); % Map of facies deposited by transport for current iteration +transFacies = uint16(0); +flowCount = uint16(0); +flowRoute = zeros(glob.ySize, glob.xSize); +transLength = uint16(0); +destCells = uint16(0); +noRoute = uint16(0); +foundDeeperCell = false; +dummy1 = uint16(0); +dummy2 = uint16(0); + +topog = glob.strata(:,:,iteration); + + + +% Calculate transport from each point on the grid +yArray=1:glob.ySize; +xArray=1:glob.xSize; +yArrayRand = yArray(randperm(length(yArray))); +xArrayRand = xArray(randperm(length(xArray))); + + + + +for i=1:glob.ySize + for j=1:glob.xSize + + y=yArrayRand(i); + x=xArrayRand(j); + + + flowThick = glob.initialTransportThickMap(y,x); % Trans thick map contains the transportable + %thicknesses calculated in the production routines, derived from the glob.erosionPerc value + prodFacies = glob.facies{y,x,iteration}(1); % Get the in-situ produced facies for this cell + + %Calculate the gradient in the specified direction + [nextX, nextY,foundDeeperCell] = currentDirection(x, y, glob,topog); + flowGradient = topog(y,x) - topog(nextY,nextX); + + % If... + % source cell contains primary producer facies 1,2, or 3 and are below sea level + % flow thickness exceeds set threshold + % flow gradient is less than the threshold + % Call the cross platform transportation + + if foundDeeperCell==true && ... + prodFacies <= glob.maxProdFacies && prodFacies > 0 && glob.wd(y,x,iteration) >= 0 && flowThick > glob.faciesThicknessPlotCutoff &&... + flowGradient <= glob.transportGradient(prodFacies) + + transFacies = glob.transportProductFacies(prodFacies); % Get the transport facies produced from the in-situ facies + + + % Remove the transportable thickness from the relevant array records + glob.strata(y,x,iteration) = glob.strata(y,x,iteration) - flowThick; + glob.thickness{y,x,iteration}(1) = glob.thickness{y,x,iteration}(1) - flowThick; + + [topog, glob] = calcOneCrossPlat(glob, iteration, topog, y, x, flowThick, prodFacies, transFacies, noRoute); + + + flowCount = flowCount + 1; + + + % If... + % source cell contains primary producer facies 1,2, or 3 and are below sea level + % flow thickness exceeds set threshold + % flow gradient is greater than the threshold + % Call the steepest descent + elseif foundDeeperCell==true && ... + prodFacies <= glob.maxProdFacies && prodFacies > 0 && glob.wd(y,x,iteration) >= 0 && flowThick > glob.faciesThicknessPlotCutoff &&... + flowGradient > glob.transportGradient(prodFacies) + + transFacies = glob.transportProductFacies(prodFacies); % Get the transport facies produced from the in-situ facies + + + % Remove the transportable thickness from the relevant array records + glob.strata(y,x,iteration) = glob.strata(y,x,iteration) - flowThick; + glob.thickness{y,x,iteration}(1) = glob.thickness{y,x,iteration}(1) - flowThick; + + [topog,glob] = calcOneTransportEvent(glob, iteration, topog, y, x, flowThick, prodFacies, transFacies, noRoute); + flowCount = flowCount + 1; + end + end +end + + +totalTransported = 0.0; + + +for y=1:glob.ySize + for x=1:glob.xSize + checkProcess=strcmp(glob.siliciclasticsRoutine,'on'); + if checkProcess==1; + %calculate the amount of siliciclastic + kl=find(glob.facies{y,x,iteration}(:)==8); + if isempty(kl) + clast = 0; + else + clast = sum(glob.thickness{y,x,iteration}(kl)); + end + else + clast =0; + end + % Increase deposited thickness by amount of deposition in all the transported + % facies at y,x + %deposited thickness is the amount of deposited without the + %siliciclastic + if glob.facies{y,x,iteration}(1)<=glob.maxProdFacies + oneThickness=sum(glob.thickness{y,x,iteration}(:))-sum(glob.thickness{y,x,iteration}(1))-clast; + else + oneThickness=sum(glob.thickness{y,x,iteration}(:))-clast; + end + + %Update the strata array + glob.strata(y,x,iteration) = glob.strata(y,x,iteration) + oneThickness; + + totalTransported = totalTransported + oneThickness; + + % Decrease WD by amount of deposition + glob.wd(y,x,iteration) = glob.wd(y,x,iteration) - oneThickness; + + if oneThickness>0; destCells=destCells+1; end + end +end + +fprintf('Trans %3.2f in %d flows to %d cells. Failed %d.', totalTransported, flowCount, destCells, noRoute); + +end + + +function [topog,glob] = calcOneCrossPlat(glob, iteration, topog, startY,startX,flowThick, prodFacies, transFacies, noRoute) +% calculates the transport from cell startX startY and records the deposited thickness. The cross platform moves material in areas with +% gradients less than the threshold. When the threshold is exceeded the cross platform terminates, the material is passed to steepeset descent +% and the cross platform cannot start again. + +% Initialize variables +flowDone = false; +flowRoute = zeros(glob.ySize, glob.xSize); +flowLength = uint16(0); % length of flow in number of grid cells +numFaciesInCell = uint8(0); +flowRoute(startY, startX) = flowLength; +steepest=false; %if steepest=false the steepest descent has not called +df=0.5; %deposited fraction + +y = startY; +x = startX; + + +%If the gradient is greater than the threshold, the steepest +%descent algorithm is called. When the steepest has been called, the +%cross platfrom algorithm cannot run again. + +while flowDone == false + + if steepest == false + + %Define current direction + %Check the next cell in the current direction and calculate the flowGradient. + + [nextX, nextY,foundDeeperCell] = currentDirection(x, y, glob,topog); + + + flowLength = flowLength + 1; + flowRoute(nextY, nextX) = flowLength; + + + %Accumulation at the production cell might create a steep gradient. + %Skip the current cell and check the gradient for the cells away from the production cell + + if y==startY && x==startX + flowGradient = 0; + else + flowGradient = topog(y,x) - topog(nextY,nextX); + end + + %If flow gradient + if flowGradient < glob.transportGradient(prodFacies) + + + %there is a deeper cell (and it is below SL) + availableSpace=glob.SL(iteration)-topog(nextY,nextX); + if availableSpace>=0 %&& foundDeeperCell==true + + if foundDeeperCell==0 + y=nextY+1; + x = nextX; + else + x = nextX; + y = nextY; + end + + if y==glob.ySize %The flow terminates when the edge of the platform has been reached. + flowDone = true; + end + + else %if the next cell is above SL or the current cell is the deepest cell. Then + %deposit everything at the previous cell, check for deposition above SL and exit. + %starting with the current cell and all tha adjacent cells + xInc = [0 1 -1 0 1 -1 0 -1 1]; + yInc = [0 0 0 -1 -1 -1 1 1 1] ; + a=1; + while flowThick > 0 %Conservation of mass. All available material must be distributed somewhere. + wrapx=x+xInc(a); + wrapy=y+yInc(a); + + %The flow terminates when the edge of the platform has benn reached at y direction. + if wrapy>glob.ySize;wrapy=glob.ySize;flowThick=0;end + if wrapy<1;wrapy=1;flowThick=0;end + + %Wrapping boundary at x direction + if wrapx>glob.xSize;wrapx=1;end + if wrapx<1;wrapx=glob.xSize;end + + %No deposition above SL + availableSpace=glob.SL(iteration)-topog(wrapy,wrapx); + if availableSpace>0 + if flowThick > availableSpace + depositHere = availableSpace; + else + depositHere = flowThick; + end + + flowThick = flowThick - depositHere; + + glob.numberOfLayers(wrapy,wrapx,iteration)=glob.numberOfLayers(wrapy,wrapx,iteration)+1; + glob.facies{wrapy, wrapx,iteration}(glob.numberOfLayers(wrapy,wrapx,iteration)) = transFacies; + glob.thickness{wrapy, wrapx,iteration}(glob.numberOfLayers(wrapy,wrapx,iteration)) = depositHere; + topog(wrapy,wrapx) = topog(wrapy,wrapx) + depositHere; % NB this means flows will interact with deposits of previous flow + + a=a+1; + if a>9 + y=y+1; + a=1; + if y>=glob.ySize + flowThick=0; + end + end + else + a=a+1; + if a>9 + y=y+1; + a=1; + if y>=glob.ySize + flowThick=0; + end + end + end + % end + % end + end + flowDone=true; + end + else % a local gradient greater than the threshold has been found. Terminate cross platform and call steepest descent. + steepest = true; + end + + else % local Gradient greater than the threshold. Call steepest. + + [topog,glob] = calcOneTransportEvent(glob, iteration, topog, y, x, flowThick, prodFacies, transFacies, noRoute); + flowDone = true; + + end +end +end + +function [topog, glob] = calcOneTransportEvent(glob, iteration, topog, startY, startX, flowThick, prodFacies, transFacies, noRoute) +% calculates the transport from cell startX startY and records the deposited thickness in +% oneTransThickMap + +% Initialize variables +flowDone = false; +flowRoute = zeros(glob.ySize, glob.xSize); +flowLength = uint16(0); % length of flow in number of grid cells +numFaciesInCell = uint8(0); +flowRoute(startY, startX) = flowLength; +flowGradient = 0; +distLimit=3; + +y = startY; +x = startX; + +while flowDone == false + + [foundDeeperCell, deepestX, deepestY] = findDeepestNeighbourCell(topog, x, y,glob); + flowGradient = topog(y,x) - topog(deepestY, deepestX); + + % Lowers cells found adjacent to the current flow cell, and gradient > minimum threshold, so deposit and + % prepare to carry on flow + + %TRANSPORT + if foundDeeperCell == true + if flowGradient > glob.transportGradient(prodFacies) + flowLength = flowLength + 1; + flowRoute(deepestY, deepestX) = flowLength; + x = deepestX; + y = deepestY; + flowDone = false; % Because there may be yet deeper cells adjacent to this one... + else + % Calculate the thickness to deposit based on a proportional decay of + % thickness and a minimum cutoff + + if flowThick < 0.01 % gradient too small and flow too thin + depositHere = flowThick; + flowThick = 0; + flowDone = true; + + else %gradient too small and thick flow + + depositHere = flowThick * 0.5; + %------- + availableSpace=glob.SL(iteration)-topog(y,x); + if availableSpace==0 + x = deepestX; + y = deepestY; + end + if depositHere>availableSpace + depositHere=availableSpace; + end + %-------- + flowThick = flowThick - depositHere; + end + + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = transFacies; + glob.thickness{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = depositHere; + topog(y, x) = topog(y, x) + depositHere; + + end + + % DEPOSIT + % Else no lower cell found, so deposit all flow at the last deepest cell x y unless it is the start cell + else if not(x == startX && y == startY) + availableSpace=glob.SL(iteration)-topog(y,x); + if flowThick<=availableSpace; + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = transFacies; + glob.thickness{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = flowThick; + topog(y, x) = topog(y, x) + flowThick; % NB this means flows will interact with deposits of previous flow + + flowDone = true; % End the flow because no deeper cell found or below gradient threshold so deposit in this low + else + %NB add this on 26/5/2015 to avoid building over sea level - estani + %& george + depositHere=availableSpace; + flowThick=flowThick-depositHere; + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = transFacies; + glob.thickness{y, x,iteration}(glob.numberOfLayers(y,x,iteration)) = depositHere; + topog(y, x) = topog(y, x) + depositHere; + + + xInc = [0 -1 1 -1 1 -1 1 0 ]; + yInc = [-1 1 -1 -1 0 0 1 1 ]; + a=1; + checkSurr=0; + + while flowThick>0; + ny=y+yInc(a) ; + if ny<1 ; ny=1; end + if ny>glob.ySize; ny=glob.ySize; end + nx=x+xInc(a); + if nx<1 ; nx=1; end + if nx>glob.xSize; nx=glob.xSize; end + availableSpace=glob.SL(iteration)-topog(ny,nx); + if availableSpace>flowThick + depositHere=flowThick; + else + depositHere=availableSpace; + end + flowThick=flowThick-depositHere; + + + if depositHere>0; + glob.numberOfLayers(ny,nx,iteration)=glob.numberOfLayers(ny,nx,iteration)+1; + glob.facies{ny,nx,iteration}(glob.numberOfLayers(ny,nx,iteration)) = transFacies; + glob.thickness{ny,nx,iteration}(glob.numberOfLayers(ny,nx,iteration)) = depositHere; + topog(ny,nx) = topog(ny,nx) + depositHere; + + end + + if availableSpace<=0; + checkSurr=1; + end + + + a=a+1; + if a>8 + if checkSurr==1 + flowThick = 0; + else + y=y+1; + a=1; + if y==glob.ySize + flowThick=0; + end + end + end + + end + flowDone=true; + end + else + flowDone = true; % Should only happen if deepest cell is start cell */ + noRoute = noRoute + 1; + end + end +end + +% fprintf('\n'); +% for k=1:glob.ySize +% for m=1:glob.xSize +% if flowRoute(k,m)>0 fprintf('(%d %d) %d %5.4f ',m,k,flowRoute(k,m), topog(k,m)); end +% end +% end +% fprintf('\n'); +end + + + +function [nextX, nextY,foundDeeperCell] = currentDirection(x, y, glob,topog) +foundDeeperCell = false; +deepestCellHeight = topog(y,x); +nextY = y; % These will be the returned values if no deeper cell is found. Important because used in gradient calc above +nextX = x; +% xInc = [0 -1 1 -1 1 -1 1 0 ]; +% yInc = [-1 1 -1 -1 0 0 1 1 ]; + xInc = [0 -1 1]; + yInc = [1 1 1 ]; + + +%dummy=0;% check if out of boundaries is the deepest cell +for k = 1:3 % Loop through the adjacent cells + + % Calculate xco and yco from the value of x y modifed by the increments in + % xInc and yInc. Should give a coord for each of the adjacent 8 cells + xWrap = x + xInc(k); + yWrap = y + yInc(k); + + checkProcess=strcmp(glob.wrapRoutine,'unwrap'); + if checkProcess==1; + %if in the edge, close boundary + if xWrap < 1; xWrap=1; end; + if xWrap > glob.xSize; xWrap=glob.xSize; end; + if yWrap < 1; yWrap = 1; end; + if yWrap > glob.ySize; yWrap = glob.ySize; end; + else + %wrap + if xWrap < 1; xWrap=glob.xSize; end; + if xWrap > glob.xSize; xWrap=1; end; + if yWrap < 1; yWrap = glob.ySize; end; + if yWrap > glob.ySize; yWrap = 1; end; + end + + % Check the current cell in the neighbour loop and see if it is the + % deepest found so far and not already visited + %Cell is deeper + if topog(yWrap, xWrap) < deepestCellHeight %&& flowRoute(yWrap,xWrap) == 0 + foundDeeperCell = true; + deepestCellHeight = topog(yWrap, xWrap); % added 5 Jan 2014 - elimnates depositional ridges bug on the slope + nextX = xWrap; + nextY = yWrap; + %two (or more cells) have the same depth. Move to the closest + %one + elseif topog(yWrap, xWrap) == deepestCellHeight + distDeep = sqrt(double((y-nextY)^2 + (x-nextX)^2)); + distNew = sqrt(double((y-yWrap)^2 + (x-xWrap)^2)); + if distNew <= distDeep + foundDeeperCell = true; + deepestCellHeight = topog(yWrap, xWrap); % added 5 Jan 2014 - elimnates depositional ridges bug on the slope + nextX = xWrap; + nextY = yWrap; + end + end +end +end + +function [foundDeeperCell, deepestX, deepestY] = findDeepestNeighbourCell(topog,x,y,glob) +n=0; +foundDeeperCell = false; +deepestCellHeight = topog(y,x); +deepestY = y; % These will be the returned values if no deeper cell is found. Important because used in gradient calc above +deepestX = x; + xInc = [0 -1 1 -1 1 -1 1 0 ]; + yInc = [-1 1 -1 -1 0 0 1 1 ]; + + +%dummy=0;% check if out of boundaries is the deepest cell +for k = 1:8 % Loop through the adjacent cells + + % Calculate xco and yco from the value of x y modifed by the increments in + % xInc and yInc. Should give a coord for each of the adjacent 8 cells + xWrap = x + xInc(k); + yWrap = y + yInc(k); + + checkProcess=strcmp(glob.wrapRoutine,'unwrap'); + if checkProcess==1; + %if in the edge, close boundary + if xWrap < 1; xWrap=1; end; + if xWrap > glob.xSize; xWrap=glob.xSize; end; + if yWrap < 1; yWrap = 1; end; + if yWrap > glob.ySize; yWrap = glob.ySize; end; + else + %wrap + if xWrap < 1; xWrap=glob.xSize; end; + if xWrap > glob.xSize; xWrap=1; end; + if yWrap < 1; yWrap = glob.ySize; end; + if yWrap > glob.ySize; yWrap = 1; end; + end + + % Check the current cell in the neighbour loop and see if it is the + % deepest found so far and not already visited + %Cell is deeper + if topog(yWrap, xWrap) < deepestCellHeight %&& flowRoute(yWrap,xWrap) == 0 + foundDeeperCell = true; + deepestCellHeight = topog(yWrap, xWrap); % added 5 Jan 2014 - elimnates depositional ridges bug on the slope + deepestX = xWrap; + deepestY = yWrap; + %two (or more cells) have the same depth. Move to the closest + %one + elseif topog(yWrap, xWrap) == deepestCellHeight + distDeep = sqrt(double((y-deepestY)^2 + (x-deepestX)^2)); + distNew = sqrt(double((y-yWrap)^2 + (x-xWrap)^2)); + if distNew <= distDeep + foundDeeperCell = true; + deepestCellHeight = topog(yWrap, xWrap); % added 5 Jan 2014 - elimnates depositional ridges bug on the slope + deepestX = xWrap; + deepestY = yWrap; + end + end +end +n=n+1; +end + + diff --git a/calculateTransportFraction.m b/calculateTransportFraction.m new file mode 100644 index 0000000..c35dea6 --- /dev/null +++ b/calculateTransportFraction.m @@ -0,0 +1,31 @@ +function [transportFraction] = calculateTransportFraction(y,x,oneFacies,glob,iteration) + +waveBase=60; %The wave base depth. Max depth where transportation occurs + +%The minimum depth for transportation. Above this depth transportation decreases to zero. +minD=zeros(4,1); + +for i=1:4 + if i==1 + minD(i)=0.5; + elseif i==2 + minD(i)=0.3; + elseif i==3 + minD(i)=0.0; + else + minD(i)=0.5; + end +end + +% %The transported fraction is zero at the surface, increases down to a min depth and decreases with depth +%and below wave base (or a user defined depth) the transported fraction becomes zero + +if glob.wd(y,x,iteration) wave break depth + k = glob.ySize; % Assumed that the waves originate at y = glob.ySize end of the grid + while k > 1 + oneFetchLength = 0; + while k >= 1 && glob.wd(k,x,iteration) > waveAmplitude(k,x) + fetchDistance(k,x) = oneFetchLength * glob.dx; + oneFetchLength = oneFetchLength + 1; + k = k - 1; + end + k = k - 1; + end + end + + % Division of zero water depth by ^-2 creates inf amplitude values so replace these with zero + waveAmplitude(isinf(waveAmplitude)) = 0; + + % Use the maximum values in each map to normalise the values + maxFetchDistance = max(max(fetchDistance)); + maxWaveAmplitude = max(max(waveAmplitude)); + normFetchDistance = fetchDistance ./ maxFetchDistance; + normWaveAmplitude = waveAmplitude ./ maxWaveAmplitude; + + % Calcluate a wave energy approximation that is a sum of the fetch and amplitude effects + glob.waveEnergy = normFetchDistance; + glob.waveEnergy(isnan(glob.waveEnergy)) = 0; % Set any NaN values to zero + + % Finally, smooth the wave energy array to represent time-averaged effects of various wave sizes, overrunning platform margin, etc +% kernel = 0.125 * ones(3); % Create a 3x3 low-pass filter kernel +% paddedWaveEnergy = zeros(glob.ySize, glob.xSize * 2); % Convolution has unwanted edge effects so need to add padding around wave energy matrix edges +% padWidth = round(glob.xSize / 2); % Define the size of the padded edges +% paddedWaveEnergy(:,1:padWidth) = repmat(glob.waveEnergy(:,1), 1, padWidth); % Add the pad on the x=1 side +% paddedWaveEnergy(:,(round(glob.xSize * 2) - padWidth + 1): round(glob.xSize*2)) = repmat(glob.waveEnergy(:,glob.xSize), 1, padWidth); % Add the pad on the x=xSize side +% paddedWaveEnergy(:,(padWidth + 1):((round(glob.xSize * 2) - padWidth))) = glob.waveEnergy; +% +% paddedWaveEnergy = conv2(paddedWaveEnergy, kernel, 'same'); +% glob.waveEnergy(1:glob.ySize, 1:glob.xSize) = paddedWaveEnergy(1:glob.ySize, (padWidth+1): (round(glob.xSize * 2) - padWidth)); + glob.waveEnergy = glob.waveEnergy ./ max(max(glob.waveEnergy)); % Renormalise 0-1 after smoothing +end \ No newline at end of file diff --git a/carboCATGUI.m b/carboCATGUI.m new file mode 100644 index 0000000..70746cf --- /dev/null +++ b/carboCATGUI.m @@ -0,0 +1,44 @@ + +% CarboCATLite is a reduced functionality version of CarboCAT2018 based on a +% version of the code from 2016 produced by Burgess, Koslowski and +% Antonatos, tidied up and recoded in parts by Burgess 2020, to run quickly +% and with a focus on 2d section output + +graph.main = 0; +graph.f1 = 0; +graph.f2 = 0; +graph.f3 = 0; +graph.f4 = 0; +graph.f5 = 0; +graph.f6 = 0; +%graph.mapMovie = zeros(1,glob.maxIts); + +% Create initial members of global structures used throughout the code, +% only those that are needed for the initial GUI display +glob.modelName = ''; +glob.xSize = 10; +glob.ySize = 100; +stats.totalFaciesVolume = 0; + +%process parameters +% glob.productionProfile='' ; +% glob.CARoutine=''; +% glob.transportationRoutine=''; +% glob.siliciclasticsRoutine=''; +% glob.concentrationRoutine=''; +% glob.soilRoutine=''; +% glob.transportFaciesRoutine=''; +% glob.refloodingRoutine=''; +% glob.seaLevelRoutine=''; +% glob.wrapRoutine=''; + +set(0,'RecursionLimit',1000); + +glob = initializeGUI(glob, stats, graph); + + + + + + + diff --git a/colorMaps/carboCAT2FaciesColours.txt b/colorMaps/carboCAT2FaciesColours.txt new file mode 100644 index 0000000..4fcda7b --- /dev/null +++ b/colorMaps/carboCAT2FaciesColours.txt @@ -0,0 +1,9 @@ + 9.0980392e-01 5.5686275e-01 5.5686275e-01 + 6.4705882e-01 8.2352941e-01 9.1372549e-01 + 5.5686275e-01 5.8823529e-02 2.0784314e-01 + 5.8823529e-02 4.0784314e-01 5.5686275e-01 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 0.0000000e+00 0.0000000e+00 diff --git a/colorMaps/carboCAT3FaciesColours.txt b/colorMaps/carboCAT3FaciesColours.txt new file mode 100644 index 0000000..111d48f --- /dev/null +++ b/colorMaps/carboCAT3FaciesColours.txt @@ -0,0 +1,9 @@ + 9.0980392e-01 5.5686275e-01 5.5686275e-01 + 6.4705882e-01 8.2352941e-01 9.1372549e-01 + 6.0000000e-01 1.0000000e+00 8.0000000e-01 + 5.5686275e-01 5.8823529e-02 2.0784314e-01 + 5.8823529e-02 4.0784314e-01 5.5686275e-01 + 0.0000000e+00 6.0000000e-01 2.9803922e-01 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 0.0000000e+00 0.0000000e+00 diff --git a/colorMaps/colorMapCA2Facies.mat b/colorMaps/colorMapCA2Facies.mat new file mode 100644 index 0000000..ffd8c19 Binary files /dev/null and b/colorMaps/colorMapCA2Facies.mat differ diff --git a/diffuseConcentration.m b/diffuseConcentration.m new file mode 100644 index 0000000..4fbd6a2 --- /dev/null +++ b/diffuseConcentration.m @@ -0,0 +1,204 @@ +function [glob]=diffuseConcentration(glob,iteration) +%diffuse + +startY=glob.yInitConcentration; +startX=glob.xInitConcentration; +diffusionOnYplus=glob.diffYplus; %of concentration, in m2/yr +diffusionOnYminus=glob.diffYminus; %of concentration, in m2/yr +diffusionOnXplus=glob.diffXplus; %of concentration, in m2/yr +diffusionOnXminus=glob.diffXminus; %of concentration, in m2/yr + +totalYears=glob.deltaT.*1000000; + +totalLoops=round(sqrt((double(glob.ySize))^2+(double(glob.xSize))^2)); +loopT=glob.deltaT*1000000/totalLoops; %in years + +%check if this values is small enough, if not decrease by an order of 10 and re-check +testDiffusion=((glob.dx)^2)/(3*loopT); + +diffArray=[diffusionOnYplus,diffusionOnYminus,diffusionOnXplus,diffusionOnXminus]; +sortedDiffArray=sort(diffArray,'descend'); + +while sortedDiffArray(1)>=testDiffusion %only check for the larger diffusion value + loopT=loopT/10; + testDiffusion=((glob.dx)^2)/(3*loopT); +end + +%get the water depth +waterDepthMap=glob.wd(:,:,iteration); +waterDepthMap(waterDepthMap<=0)=0; + +%calculate gradient + + +%loop through cells around radius until cell underwater is found, in case +%the source point has been filled + +foundCell=false; + +r=0; +leaveOut=0; + + +while foundCell==false + + + [yArray,xArray,length]=calculateNeighbourCellsFromRadius(r,leaveOut); + + for yPoint=1:length + + for xPoint=1:length + + yco=startY+yArray(yPoint); + xco=startX+xArray(xPoint); + + if yco<1; yco=1; end + if yco>glob.ySize; yco=glob.ySize; end + + if xco<1; xco=1; end + if xco>glob.xSize; xco=glob.xSize; end + + %check if cell is underwater + oneWaterDepth=waterDepthMap(yco,xco); + if oneWaterDepth>0 + foundCell=true; + newY=yco; + newX=xco; + end + + end + end + + r=r+1; + leaveOut=leaveOut+1; + +end + + +%convert Rate to Volume +glob.volumeCarbonate=glob.inputRateCarbonate.*glob.deltaT.*1000000; +%convert the map from the previous itereation to concentration +initialConcentration=glob.carbonateVolMap./waterDepthMap; +initialConcentration(waterDepthMap<=0)=0; +%calculate the concentration in the initial point +%glob.inputConcentration=glob.volumeCarbonate/waterDepthMap(newY,newX); + +totalLoops=totalYears/loopT; + +%add input concentration to the grid +%initialConcentration(newY,newX)=glob.inputConcentration+initialConcentration(newY,newX); +inputVolAtEachLoop=glob.volumeCarbonate./totalLoops; + + +%re-scale de diffusion coeff to time step + +diffusionOnYplusScaled=diffusionOnYplus.*loopT; +diffusionOnYminusScaled=diffusionOnYminus.*loopT; +diffusionOnXplusScaled=diffusionOnXplus.*loopT; +diffusionOnXminusScaled=diffusionOnXminus.*loopT; + +%totalLoops=50; +for iter=1:totalLoops + %figure(12); surf(initialConcentration); view([0 90]); caxis([0 glob.inputConcentration]); + + %add input concentration to the grid + initialConcentration(newY,newX)=inputVolAtEachLoop/waterDepthMap(newY,newX)+initialConcentration(newY,newX); + + concxminus= circshift(initialConcentration,[0,-1]); + concxplus= circshift(initialConcentration,[0,1]); + concyminus= circshift(initialConcentration,[-1,0]); + concyplus= circshift(initialConcentration,[1,0]); + + wdxminus= circshift(waterDepthMap,[0,-1]); + wdxplus= circshift(waterDepthMap,[0,1]); + wdyminus= circshift(waterDepthMap,[-1,0]); + wdyplus= circshift(waterDepthMap,[1,0]); + + nxminus=zeros(glob.ySize,glob.xSize); + nxplus=zeros(glob.ySize,glob.xSize); + + nyminus=zeros(glob.ySize,glob.xSize); + nyplus=zeros(glob.ySize,glob.xSize); + + concxminus(waterDepthMap<=0)=0; + concxplus(waterDepthMap<=0)=0; + + concyminus(waterDepthMap<=0)=0; + concyplus(waterDepthMap<=0)=0; + + %control the diffusion towards exposed cells + nxminus(wdxminus>0)=nxminus(wdxminus>0)+1; + nxplus(wdxplus>0)=nxplus(wdxplus>0)+1; + + nyminus(wdyminus>0)=nyminus(wdyminus>0)+1; + nyplus(wdyplus>0)=nyplus(wdyplus>0)+1; + + + dcdx2minus = concxminus - (nxminus.*initialConcentration); + dcdx2plus = concxplus - (nxplus.*initialConcentration); + + dcdy2minus = concyminus - (nyminus.*initialConcentration); + dcdy2plus = concyplus - (nyplus.*initialConcentration); + + %diffusion acts slower when depth increases, and faster when depth + %decreases + diffusionOnXminusScaledMap=diffusionOnXminusScaled./waterDepthMap; + diffusionOnXminusScaledMap(waterDepthMap<1.0)=diffusionOnXminusScaled; + diffusionOnXplusScaledMap=diffusionOnXplusScaled./waterDepthMap; + diffusionOnXplusScaledMap(waterDepthMap<1.0)=diffusionOnXminusScaled; + diffusionOnYminusScaledMap=diffusionOnYminusScaled./waterDepthMap; + diffusionOnYminusScaledMap(waterDepthMap<1.0)=diffusionOnXminusScaled; + diffusionOnYplusScaledMap=diffusionOnYplusScaled./waterDepthMap; + diffusionOnYplusScaledMap(waterDepthMap<1.0)=diffusionOnXminusScaled; + + deltaConcentrationxminus = (diffusionOnXminusScaledMap.*dcdx2minus) ./ ((2*glob.dx)^2); + deltaConcentrationxplus = (diffusionOnXplusScaledMap.*dcdx2plus) ./ ((2*glob.dx)^2); + + + deltaConcentrationyminus = (diffusionOnYminusScaledMap.*dcdy2minus) ./ ((2*glob.dx)^2); + deltaConcentrationyplus = (diffusionOnYplusScaledMap.*dcdy2plus) ./ ((2*glob.dx)^2); + + %re-scale the concentration to the new depth + %get the volume that has been transferred + deltaConcentrationxminus(deltaConcentrationxminus>0) = deltaConcentrationxminus(deltaConcentrationxminus>0).*wdxminus(deltaConcentrationxminus>0)./waterDepthMap(deltaConcentrationxminus>0); + deltaConcentrationxplus(deltaConcentrationxplus>0) = deltaConcentrationxplus(deltaConcentrationxplus>0).*wdxplus(deltaConcentrationxplus>0)./waterDepthMap(deltaConcentrationxplus>0); + + + deltaConcentrationyminus(deltaConcentrationyminus>0) = deltaConcentrationyminus(deltaConcentrationyminus>0).*wdyminus(deltaConcentrationyminus>0)./waterDepthMap(deltaConcentrationyminus>0); + deltaConcentrationyplus(deltaConcentrationyplus>0) = deltaConcentrationyplus(deltaConcentrationyplus>0).*wdyplus(deltaConcentrationyplus>0)./waterDepthMap(deltaConcentrationyplus>0); + + + deltaConcentration = deltaConcentrationyminus+deltaConcentrationxminus+deltaConcentrationyplus+deltaConcentrationxplus; + initialConcentration = initialConcentration+deltaConcentration; +end + +%transform to volume +glob.carbonateVolMap=initialConcentration.*waterDepthMap; + +%in case you want to plot the concentration +% figure(10) +% subplot(1,3,1); +% hold on +% % surf(glob.strata(:,:,1),'edgecolor', 'none'); view([0 90]); axis square, axis tight ;colorbar ; +% contourf(glob.strata(:,:,1),10); view([0 90]); axis square, axis tight ;colorbar ; +% contour(glob.strata(:,:,1),[0 0],'w','linewidth',2); +% xlabel(' distance (km)'); +% ylabel('distance (km)'); +% title('elevation (m)') +% subplot(1,3,2); +% hold on +% contourf(initialConcentration,5); view([0 90]); colorbar; axis square, axis tight; +% contour(glob.strata(:,:,1),[0 0],'w','linewidth',2); +% xlabel(' distance (km)'); +% ylabel('distance (km)'); +% title('carbonate concentration (m3/m)') +% subplot(1,3,3); +% hold on +% contourf(glob.carbonateVolMap,5); view([0 90]); colorbar; axis square, axis tight; +% contour(glob.strata(:,:,1),[0 0],'w','linewidth',2); +% xlabel(' distance (km)'); +% ylabel('distance (km)'); +% title('dissolved carbonate (m3)') +end + + diff --git a/diffuseSiliciclastic.m b/diffuseSiliciclastic.m new file mode 100644 index 0000000..231197a --- /dev/null +++ b/diffuseSiliciclastic.m @@ -0,0 +1,256 @@ +function [glob]=diffuseSiliciclastic(glob,iteration) +%diffuse + +topog = glob.strata(:,:,iteration-1); +availableSpace=glob.SL(iteration)-topog(:,:); + +glob.productionBySiliciclasticMap = zeros(glob.ySize,glob.xSize,glob.maxProdFacies); + +%The coordinates of the point source. Either directly or from the file +startY=glob.yInitSili; +startX=glob.xInitSili;% point source +sLength = glob.sourceLength; %line source + +%Concentration rate. It describes how far the siliciclastic will spread in +%each direction. High values means it goes farther. +diffusionOnY=glob.diffY; %of concentration, in m2/yr +diffusionOnX=glob.diffX; %of concentration, in m2/yr + +totalYears=glob.deltaT.*1000000; + +loopT=100; %in years, %defines the number of difussion steps for each iteration. High steps more smooth deposition + +%check if this values is small enough, if not decrease by an order of 10 +%and re-check. Finite difference stability. +testDiffusion=((glob.dx)^2)/(3*loopT); +while diffusionOnY>=testDiffusion +loopT=loopT/10; +testDiffusion=((glob.dx)^2)/(3*loopT); +end + +while diffusionOnX>=testDiffusion +loopT=loopT/10; +testDiffusion=((glob.dx)^2)/(3*loopT); +end + +%get the water depth +waterDepthMap=glob.wd(:,:,iteration); + +waterDepthMap(waterDepthMap<=0)=0; + + +%loop through cells around radius until cell underwater is found, in case +%the source point has been filled + +foundCell=false; + +r=0; +leaveOut=0; + + +while foundCell==false + + + [yArray,xArray,length]=calculateNeighbourCellsFromRadius(r,leaveOut); + + for yPoint=1:length + + for xPoint=1:length + + yco=startY+yArray(yPoint); + xco=startX+ceil(sLength/2)+xArray(xPoint); + + if yco<1; yco=1; end + if yco>glob.ySize; yco=glob.ySize; end + + if xco<1; xco=1; end + if xco>glob.xSize; xco=glob.xSize; end + + %check if cell is underwater + oneWaterDepth=glob.wd(yco,xco); + if oneWaterDepth>0 + foundCell=true; + newY=yco; +% newX=xco; + end + + end + end + + r=r+1; + leaveOut=leaveOut+1; + +end + +%Keep the coordinates of the source constant. +newX = startX; + + + +%get the initial concentration from the previous time-step or zero +initialConcentration=zeros(glob.ySize,glob.xSize); + +%add input concentration to the grid +initialConcentration(newY,newX:newX+sLength)=glob.inputSili; + + + + totalLoops=totalYears/loopT; + if totalLoops>5000;totalLoops=5000;end %not necessary. Time consuming is totalLoops>10.000 + + + %re-scale de diffusion coeff to time step + + diffusionOnYScaled=diffusionOnY.*loopT; + diffusionOnXScaled=diffusionOnX.*loopT; + + + for iter=1:totalLoops + + dcdx2 = circshift(initialConcentration,[0,-1]) + circshift(initialConcentration,[0,1]) - (2.*initialConcentration);%propagation in -x and +x direction +% dcdx2 = circshift(initialConcentration,[0,-1]) - (1.*initialConcentration);%propagation only in -x or +x direction +% dcdy2 = circshift(initialConcentration,[-1,0]) + circshift(initialConcentration,[1,0]) - (2.*initialConcentration);%propgation in -y and +y direction + dcdy2 = circshift(initialConcentration,[1,0]) -(1.*initialConcentration);%propgation only in -y or +y direction + + + + deltaConcentrationx = (diffusionOnXScaled.*dcdx2) ./ ((2*glob.dx)^2);%diffusion along x-axis + deltaConcentrationy = (diffusionOnYScaled.*dcdy2) ./ ((2*glob.dx)^2);%diffusion along y-axis + + %correct for material beeing diffused over cells above water + + for y=1:glob.ySize + for x=1:glob.xSize + + + if waterDepthMap(y,x)<=0 + + + if deltaConcentrationx(y,x)> 0 %check for concentration shifts in the x direction + + %check the surrounding cells + + yco(1)=y; + yco(2)=y; + xco(1)=x-1; if xco(1)<1; xco(1)=1; end + xco(2)=x+1; if xco(2)>glob.xSize; xco(2)=glob.xSize; end + + conc(1)=initialConcentration(yco(1),xco(1)); + conc(2)=initialConcentration(yco(2),xco(2)); + + + %put the concentration back into original cell + differenceBetweenConcentrations=initialConcentration(y,x)-conc; + diffuseBack = (-diffusionOnXScaled.*differenceBetweenConcentrations) ./ ((2*glob.dx)^2); + deltaConcentrationx(y,x)=0; + deltaConcentrationx(yco(1),xco(1))=deltaConcentrationx(yco(1),xco(1))+diffuseBack(1); + deltaConcentrationx(yco(2),xco(2))=deltaConcentrationx(yco(2),xco(2))+diffuseBack(2); + + + end + if deltaConcentrationy(y,x)>0%check for concentration shifts in the y direction + + %check the surrounding cells + + xco(1)=x; + xco(2)=x; + yco(1)=y-1; if yco(1)<1; yco(1)=1; end + yco(2)=y+1; if yco(2)>glob.ySize; yco(2)=glob.ySize; end + + conc(1)=initialConcentration(yco(1),xco(1)); + conc(2)=initialConcentration(yco(2),xco(2)); + + + %put the concentration back into original cell + differenceBetweenConcentrations=initialConcentration(y,x)-conc; + diffuseBack = (-diffusionOnYScaled.*differenceBetweenConcentrations) ./ ((2*glob.dx)^2); + + deltaConcentrationy(y,x)=0; + deltaConcentrationy(yco(1),xco(1))=deltaConcentrationy(yco(1),xco(1))+diffuseBack(1); + deltaConcentrationy(yco(2),xco(2))=deltaConcentrationy(yco(2),xco(2))+diffuseBack(2); + + + end + end + end + end + + + + %change the concentration by the deltaConcentration + deltaConcentration = deltaConcentrationy+deltaConcentrationx; + initialConcentration = initialConcentration+deltaConcentration; + + + + end + + + %delete concetrations below the facies cut off thickness + initialConcentration(initialConcentration 0 + glob.numberOfLayers(y,x,iteration)=glob.numberOfLayers(y,x,iteration)+1; + glob.facies{y,x,iteration}(glob.numberOfLayers(y,x,iteration)) = 8; + glob.thickness{y,x,iteration}(glob.numberOfLayers(y,x,iteration)) = glob.concentration(y,x,iteration); + end + + if glob.concentration(y,x,iteration) > availableSpace(y,x) + glob.strata(y,x,iteration) = glob.strata(y,x,iteration-1)+availableSpace(y,x); + glob.thickness{y,x,iteration}(glob.numberOfLayers(y,x,iteration)) = availableSpace(y,x); + else + glob.strata(y,x,iteration) = glob.strata(y,x,iteration-1) + glob.concentration(y,x,iteration); + end + + end + end + + %update the water depth array + glob.wd(:,:,iteration) = glob.wd (:,:,iteration)- glob.concentration(:,:,iteration); + glob.wd(glob.wd<0)=0; + +%calculate the effect od siliciclastic to production. +%The effect is calculated as the ratio of the value at each cell over the +%value that kills production completely (kill value) +%Higher kill value =lower effect of the small concentrations, siliciclastic effect is apparent only close to the source +%Lower kill value = higher effect of small concentrations, siliciclastic effect is apparent farther from the source + +%The effect of siliciclastic on each carbonate factory is quantified by the +%kill value for the factory. +killValue1 = 1; %siliciclastic effect on factory 1 +killValue2 = .5; %siliciclastic effect on factory 2 +killValue3 = .1; %siliciclastic effect on factory 3 + +for i=1:glob.maxProdFacies + if i==1 + glob.productionBySiliciclasticMap(:,:,i) = 1.-(initialConcentration / killValue1); + elseif i==2 + glob.productionBySiliciclasticMap(:,:,i) = 1.-(initialConcentration / killValue2); + else + glob.productionBySiliciclasticMap(:,:,i) = 1.-(initialConcentration / killValue3); + end +glob.productionBySiliciclasticMap(glob.productionBySiliciclasticMap<0)=0; + + %calculate scale from 0-1 + +% a=1/(glob.maxConcValue-glob.minConcValue); +% b=-a*glob.minConcValue; +% +% glob.productionByConcentrationMap=a.*initialConcentration+b; +% +% glob.productionByConcentrationMap(glob.productionByConcentrationMap<0)=0; +% glob.productionByConcentrationMap(glob.productionByConcentrationMap>1)=1; + +end + + diff --git a/functions list.xlsx b/functions list.xlsx new file mode 100644 index 0000000..8d1cabe Binary files /dev/null and b/functions list.xlsx differ diff --git a/getData.m b/getData.m new file mode 100644 index 0000000..c0e563e --- /dev/null +++ b/getData.m @@ -0,0 +1,176 @@ +function glob = getData(glob,stats) +%get data from the model, use this function for taking variables outside +%the GUI workspace, otherwise unavailable from the command window +%for example +%assignin('base','strata',glob.strata(:,:,1)); +%puts on the workspace the variable strata with the values from +%glob.strata(:,:,1) +nameLogs={'w11-3' 'w11-1' 'w11-2' 'Durslton Head' 'Swanworth Quarry' 'Hells Bottom (incomplete)' 'Worborrow Tout' 'Mupe Bay' 'Poxwell (incomplete)' 'North Portland' 'South Portland' 'Portesham'}; +xpoint=[10 12 13 25 37 46 55 64 81 88 90 103 0]; +ypoint=[26 29 35 41 39 35 36 36 29 49 57 24 0]; +numberOfLocations=size(xpoint,2)-1; +% +% for it=1:glob.totalIterations +% for locations=1:numberOfLocations-1 +% waterDepthATPoint(it,1,locations)=glob.wd(ypoint(locations),xpoint(locations),it); +% +% if glob.numberOfLayers(ypoint(locations),xpoint(locations),it)>0 +% for layer=1:glob.numberOfLayers(ypoint(locations),xpoint(locations),it) +% faciesATPoint(it,layer,locations)=glob.facies{ypoint(locations),xpoint(locations),it}(layer); +% thicknessATPoint(it,layer,locations)=glob.thickness{ypoint(locations),xpoint(locations),it}(layer); +% end +% else +% faciesATPoint(it,1,locations)=0; +% thicknessATPoint(it,1,locations)=0; +% end +% end +% end +% +% assignin('base','faciesATPoint',faciesATPoint); +% assignin('base','thicknessATPoint',thicknessATPoint); +% assignin('base','waterDepthATPoint',waterDepthATPoint); + + + +[ average] = calculateSurfaceVariation( glob ); +%average=average/0.007; +assignin('base','averageSurfVarA',average); + +for t=2:glob.totalIterations +% totalThickness=glob.strata(:,:,t)-glob.strata(:,:,t-1); +% wd=glob.wd(:,:,t); +% wd(wd<0)=0; +% wdprev=glob.wd(:,:,t); +% wdprev(wdprev<0)=0; +% dwd=wd; +% accomm(:,:,t) = totalThickness + wd; +% %dWaterD=dWaterD+(glob.wd(:,:,t)-glob.wd(:,:,t-1)); +% %accommodation = dWaterD+totalThickness(:,:,t); + accomm(:,:,t)=(glob.wd(:,:,t)-glob.wd(:,:,t-1))+(glob.strata(:,:,t)-glob.strata(:,:,t-1)); +end +wd=glob.wd(:,:,1); +wd(wd<0)=0; +totalAccommodation=sum(accomm(:)) + sum(wd(:)); +%Calculate thickness +initStrata=glob.strata(:,:,1); +finalStrata=glob.strata(:,:,glob.totalIterations); +thickness=sum(sum(finalStrata-initStrata)); + +thicknessIndex=thickness/totalAccommodation; +assignin('base','thicknessIndexA',thicknessIndex); + + +%calculate transitions +maxTransitions=double(glob.ySize) * double(glob.xSize) * double(glob.totalIterations); +transitions=0; +for t = 2: glob.totalIterations + for y=1:glob.ySize + for x=1:glob.xSize + prodFaciesPrev=glob.facies{y,x,t-1}(1); + if prodFaciesPrev>glob.maxProdFacies + prodFaciesPrev=0; + end + prodFacies=glob.facies{y,x,t}(1); + if prodFacies>glob.maxProdFacies + prodFacies=0; + end + if prodFacies~=prodFaciesPrev + transitions=transitions+1; + end + end + end +end + +transitionsIndex=transitions/maxTransitions; + +assignin('base','transitionsIndexA',transitionsIndex); + +%calculate heterogeneity fortwo facies +transitionsIndex=transitions/maxTransitions; + + +countF1=0; +countAll=0; +for t = 1: glob.totalIterations + for y=1:glob.ySize + for x=1:glob.xSize + prodFacies=glob.facies{y,x,t}(1); + if prodFacies==1 + %count facies 1 + countF1=countF1+1; + countAll=countAll+1; + else + if prodFacies==2 + countAll=countAll+1; + end + end + end + end +end + +fractionFacies=countF1/countAll; + +assignin('base','fractionFaciesA',fractionFacies); + +% figure (20) +% for a=1:12; +% hold on +% if a<=6; +% colorArray=[a/6 0 (6-a)/6]; +% else +% colorArray=[(a-6)/6 1 (6-(a-6))/6]; +% end +% plot(waterDepthATPoint(:,1,a),'color',colorArray); +% +% end +% legend(nameLogs); + + +%calculate the simillarity between measured and calculated logs + +% thicknessAndFraction=load('carbo-CAT/params/Purbeck/logsThickness.txt'); +% posit=[1 5 9 13]; +% errorFraction=zeros(numberOfLocations,1); +% for logs=1:numberOfLocations +% x=xpoint(logs); +% y=ypoint(logs); +% %Calculate thickness of calculated +% initStrata=glob.strata(y,x,1); +% finalStrata=glob.strata(y,x,glob.totalIterations); +% totalThickness=finalStrata-initStrata; +% %get thickness of measured +% minThick= int16.empty(4,0); +% maxThick= int16.empty(4,0); +% minFr= int16.empty(4,0); +% maxFr= int16.empty(4,0); +% for p=1:4 +% minThick(p)=thicknessAndFraction(logs,posit(p)); +% %max values +% maxThick(p)=thicknessAndFraction(logs,posit(p)+1); +% %min fr +% minFr(p)=thicknessAndFraction(logs,posit(p)+2); +% %max fr +% maxFr(p)=thicknessAndFraction(logs,posit(p)+3); +% end +% totalMaxThick=sum(maxThick(:)); +% totalMinThick=sum(minThick(:)); +% +% +% if totalThickness>totalMaxThick +% if totalMaxThick<999 +% errorFraction(logs)=(totalThickness-totalMaxThick)/totalMaxThick; +% end +% else +% if totalThickness glob.maxProdFacies + fprintf('WARNING: %d producing facies in the initial facies map from %s versus %d maximum facies in from the parameter file\n',checkMaxFacies, glob.initFaciesFName, glob.maxProdFacies); + end + + for f = 1:9 + faciesTotals(f) = sum(oneFaciesMap(:) == f); + fprintf('Facies %d total initial cells is %d\n', f, faciesTotals(f)); + end + + for y=1:glob.ySize + for x= 1:glob.xSize + glob.facies{y,x,1}(1) = oneFaciesMap(y,x); + if glob.facies{y,x,1}(1)>0 + glob.numberOfLayers(y,x,1)=1; + else + glob.numberOfLayers(y,x,1)=0; + end + end + end + + % Load the initial bathymetry map using the file name specified in the parameter file + oneBathymetryMap = (load(glob.initBathymetryFName, '-ASCII')); + mapSizeCheck = size(oneBathymetryMap); + if mapSizeCheck(1) ~= glob.ySize || mapSizeCheck(2) ~= glob.xSize + fprintf('WARNING: initial bathymetry map is %d %d in size, different to the %d %d model grid\n', mapSizeCheck(1:2), glob.ySize, glob.xSize); + end + + % Set the initial water depth and the elevation of the initial stratal surface to the initial water depth + % Note that zero is the model datum so initial surface elevation is zero - initial water depth + glob.wd(:,:,1) = oneBathymetryMap; + glob.strata(:,:,1) = -glob.wd(:,:,1); + glob.strataOriginalPosition = glob.strata; + + % Load a subsidence map from an external file. + oneSubsidenceMap = (load(glob.subsidenceFName, '-ASCII')); % The path to the file with the subsidence map + glob.subRateMap = oneSubsidenceMap; + glob.subRateMap = glob.subRateMap * glob.deltaT;% Adjust production rates for timestep + + if strcmp(glob.seaLevelRoutine,'file') ==0 + % initialize sea-level curve + j = 1:glob.maxIts; % implicit loop on iterations + emt = double(j) * glob.deltaT; % create a vector of emt values + glob.SL = ((sin(pi*((emt/glob.SLPeriod1)*2)))* glob.SLAmp1)+ (sin(pi*((emt/glob.SLPeriod2)*2)))* glob.SLAmp2; + else + % read the SL/water level curve and initialise the curve array in glob + imported = importdata(glob.SLCurveFName); + glob.SL = imported; + end + + % read the carbonate production rate time curve and initialise the curve array in glob + imported = importdata(glob.carbProdCurveFName); + glob.carbProdCurve = imported; + + % Correct the transport gradient with the cell size information + glob.transportGradient(:)=glob.transportGradient(:).*glob.dx./1000; + + % Set up the CA iteration counting arrays + glob.CADtCount=zeros(glob.ySize, glob.xSize); % Set Dt counter to zero + glob.CADtPerIteration = ones(glob.ySize, glob.xSize) * glob.CADtMax; % Set the timesteps required per iteration to the input parameter value + + % Set up the wave energy array map + glob.waveEnergy = zeros(glob.ySize, glob.xSize); + + %Initialize the carbonate concentration array + %Read the concentration data + fileInConcentration = fopen(glob.concFName); + glob.inputRateCarbonate=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.yInitConcentration=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.xInitConcentration=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.diffYplus=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.diffYminus=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.diffXplus=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.diffXminus=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.initialVolumeValue=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + + glob.carbonateVolMap=zeros(glob.ySize,glob.xSize); + %set an homogeneus initial concentration value (sea water is typical 0.15kg/m3) + glob.carbonateVolMap(:,:)=glob.initialVolumeValue; + + % this is for the diffusion silicilastic model + glob.productionBySiliciclasticMap=ones(glob.ySize,glob.xSize,glob.maxProdFacies); + %Initialize the siliciclastics array + fileInConcentration = fopen(glob.siliFName); + glob.inputSili=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + + glob.yInitSili=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.xInitSili=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + + glob.sourceLength=fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + + glob.diffY = fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + glob.diffX = fscanf(fileInConcentration,'%f', 1); + dummyLabel = fgetl(fileInConcentration); + + % Finally, load a colour map for the CA facies and hiatuses + glob.faciesColours = load(glob.faciesColoursFilename); + + % Set producing facies controls that depend on number of neighbour parameters + oneFacies = 1:glob.maxProdFacies; + glob.prodScaleMin(oneFacies) = glob.CARules(oneFacies,2); + glob.prodScaleMax(oneFacies) = glob.CARules(oneFacies,3); + glob.prodScaleOptimum(oneFacies) = ((glob.prodScaleMax(oneFacies) - glob.prodScaleMin(oneFacies)) / 2)+glob.prodScaleMin(oneFacies); + glob.prodWaveOptimum(oneFacies) = glob.prodWaveThresholdLow(oneFacies) + ((glob.prodWaveThresholdHigh(oneFacies) - glob.prodWaveThresholdLow(oneFacies)) / 2.0); + + if size(glob.faciesColours,1) < glob.maxFacies % So if too few rows in the colour map, give a warning... + fprintf('Only %d colours in colour map colorMaps/faciesColourMap.txt but %d facies in model\n', size(1), glob.maxFacies); + fprintf('This could get MESSY!\n\n\n\n'); + end + + stats.totalFaciesVolume = zeros(glob.maxIts, glob.maxFacies); +end + + + + + + diff --git a/inputOneModelParams.m b/inputOneModelParams.m new file mode 100644 index 0000000..bbd3aa1 --- /dev/null +++ b/inputOneModelParams.m @@ -0,0 +1,200 @@ +function [glob, inputSuccessFlag] = inputOneModelParams(glob, fName, fPName) + + inputSuccessFlag = 1; % Assume initialisation works, unless the flag is reset to zero below, for example by file read error + + %read the processes file + fileIn = fopen(fPName); + if (fileIn < 0) + fprintf('\n WARNING: file %s not found, code about to terminate\n', fPName); + inputSuccessFlag = 0; + else + fprintf('\n Reading parameters from filename %s\n', fPName); + end + + glob.CARoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.transportationRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.waveRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.siliciclasticsRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.concentrationRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.soilRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.refloodingRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.seaLevelRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + glob.wrapRoutine = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + + fprintf('CA routine is %s based \n', glob.CARoutine); + fprintf('Transport routine is %s \n', glob.transportationRoutine); + fprintf('Wave routine is %s \n', glob.waveRoutine); + fprintf('Siliciclastics are %s \n', glob.siliciclasticsRoutine); + fprintf('Carbonate concentration in water is %s \n', glob.concentrationRoutine); + fprintf('Soild deposition is %s \n', glob.soilRoutine); + fprintf('Reflooded cells are %s \n', glob.refloodingRoutine); + fprintf('Sea-level curve from %s \n', glob.seaLevelRoutine); + fprintf('Model edges are treated as %s \n \n', glob.wrapRoutine); + + fileIn = fopen(fName); + if (fileIn < 0) + fprintf('WARNING: file %s not found, code about to terminate\n', fName); + else + fprintf(' Reading parameters from filename %s\n', fName); + end + + glob.modelName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); % Read to the end of the line to skip any label text + + % Read parameters from the main parameter values file + glob.xSize = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.ySize = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.dx = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.totalIterations = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.deltaT = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + + glob.SLPeriod1 = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.SLAmp1 = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.SLPeriod2 = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.SLAmp2 = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + + glob.CADtMin = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.CADtMax = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + + glob.BathiLimit = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + + glob.maxProdFacies = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + + for j = 1:glob.maxProdFacies + glob.prodRate(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + + glob.surfaceLight(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + if glob.surfaceLight(j) > 500 + glob.extinctionCoeff(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.saturatingLight(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + else + glob.profCentre (j)= glob.surfaceLight(j); + glob.profWidth (j)= fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.profSlope (j)= fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + end + + fprintf('Facies %d produced at %3.2f m/My\n', j, glob.prodRate(j)); + + glob.transportProductFacies(j) = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.transportFraction(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.transportGradient(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.transContinueProb(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + fprintf('Breakdown of facies %d creates facies %d at rate %3.2f of accumulated thickness \n', j, glob.transportProductFacies(j), glob.transportFraction(j)); + + glob.prodWaveThresholdLow(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + glob.prodWaveThresholdHigh(j) = fscanf(fileIn,'%f', 1); + dummyLabel = fgetl(fileIn); + fprintf('Facies %d produced from wave energy %4.3f to %4.3f\n', j, glob.prodWaveThresholdLow(j), glob.prodWaveThresholdHigh(j)); + + glob.prodRate(j) = glob.prodRate(j) * glob.deltaT; % Adjust production rates for timestep + + % Calculate the water depth cutoff below which production rate is effectively zero + % Factory types will only occur above this water depth cutoff + wd = 0.0; + if glob.surfaceLight(j)>500 + while tanh((glob.surfaceLight(j) * exp(-glob.extinctionCoeff(j) * wd))/ glob.saturatingLight(j)) > 0.000001 && wd < 10000 + glob.prodRateWDCutOff(j) = wd; + wd = wd + 0.1; + end + else + while (1/ (1+((wd-glob.profCentre(j))./glob.profWidth(j)).^(2.*glob.profSlope(j))) ) > 0.000001 && wd < 10000 + glob.prodRateWDCutOff(j) = wd; + wd = wd + 0.1; + end + end + fprintf('Facies %d has production cutoff at %3.2f m water depth\n', j, glob.prodRateWDCutOff(j)); + end + + glob.subsidenceFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Subsdience map filename %s \n', glob.subsidenceFName); + + glob.CARulesFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('CA rules filename %s \n', glob.CARulesFName); + + glob.initFaciesFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Initial condition facies map filename %s \n', glob.initFaciesFName); + + glob.initBathymetryFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Initial bathymetry map filename %s \n', glob.initBathymetryFName); + + glob.concFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Carbonate concentration filename %s \n', glob.concFName); + + glob.siliFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Carbonate siliciclastic supply filename %s \n', glob.siliFName); + + glob.SLCurveFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Sea-level curve filename %s \n', glob.SLCurveFName); + + glob.carbProdCurveFName = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Maximum carbonate production rate time curve filename %s \n', glob.carbProdCurveFName); + + glob.faciesColoursFilename = fscanf(fileIn,'%s', 1); + dummyLabel = fgetl(fileIn); + fprintf('Model facies colour map filename %s \n', glob.faciesColoursFilename); + + % Read the cellular automata rules from file name in glob.CARulesFName + if exist(glob.CARulesFName, 'file') == 2 + import = importdata(glob.CARulesFName,' ',1); + glob.CARules = import.data; + else + message = sprintf('Could not find %s\nCheck path and filename exist.', glob.CARulesFName); + h1 = msgbox(message,'Loading CA rules not successful'); + inputSuccessFlag = 0; + end + + % Read the number and ages of time lines to be plotted on cross sections. Age = iteration number + glob.timeLineCount = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.timeLineAge = zeros(1,glob.timeLineCount+1); + + glob.timeLineAge = fscanf(fileIn,'%d', glob.timeLineCount); % reads glob.timeLineCount values from the file + dummyLabel = fgetl(fileIn); + fprintf('Plotting %d timelines from iteration %d to %d\n', glob.timeLineCount, glob.timeLineAge(1), glob.timeLineAge(glob.timeLineCount)); + + % Finally, read the number and ages of maps to be plotted in the relevant figure. Age = iteration number + glob.mapCount = fscanf(fileIn,'%d', 1); + dummyLabel = fgetl(fileIn); + glob.mapAge = zeros(1,glob.mapCount+1); +end \ No newline at end of file diff --git a/params/DbPlatform/CARules3_4_10_6_10.txt b/params/DbPlatform/CARules3_4_10_6_10.txt new file mode 100644 index 0000000..4d3e1b1 --- /dev/null +++ b/params/DbPlatform/CARules3_4_10_6_10.txt @@ -0,0 +1,4 @@ +Dist MinNeighbs MaxNeighbs MinTrigger MaxTrigger TriggerFacies +2 2 10 4 10 1 +2 4 10 6 10 2 +2 4 10 6 10 3 diff --git a/params/DbPlatform/carbProdConst1000Iterations.txt b/params/DbPlatform/carbProdConst1000Iterations.txt new file mode 100644 index 0000000..59fbfc9 --- /dev/null +++ b/params/DbPlatform/carbProdConst1000Iterations.txt @@ -0,0 +1,1000 @@ +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 diff --git a/params/DbPlatform/initialFaciesMapRandomY150X10.txt b/params/DbPlatform/initialFaciesMapRandomY150X10.txt new file mode 100644 index 0000000..1e4d441 --- /dev/null +++ b/params/DbPlatform/initialFaciesMapRandomY150X10.txt @@ -0,0 +1,150 @@ + 0.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 + 0.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 + 1.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 + 3.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 + 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 + 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 + 1.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 3.0000000e+00 + 1.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 + 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 + 1.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 + 3.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 + 1.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 + 0.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 + 0.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 + 0.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 + 0.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 + 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 + 1.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 + 1.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 + 0.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 + 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 + 1.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 + 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 + 0.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 + 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 + 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 + 1.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 3.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 + 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 + 1.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 + 3.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 + 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 + 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 + 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 + 3.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 + 3.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 3.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 + 3.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 + 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 + 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 + 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 + 1.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 + 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 + 0.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 + 2.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 + 1.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 0.0000000e+00 + 3.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 + 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 0.0000000e+00 3.0000000e+00 3.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 + 0.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 + 3.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 + 3.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 + 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 + 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 + 3.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 + 0.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 + 3.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 + 2.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 + 3.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 + 1.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 + 2.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 + 2.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 + 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 + 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 + 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 1.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 + 0.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 + 3.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 + 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 3.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 + 2.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 + 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 3.0000000e+00 0.0000000e+00 + 0.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 + 1.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 + 2.0000000e+00 3.0000000e+00 3.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 3.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 2.0000000e+00 0.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 0.0000000e+00 0.0000000e+00 3.0000000e+00 + 3.0000000e+00 0.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 + 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 + 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 + 3.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 3.0000000e+00 0.0000000e+00 + 1.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 + 0.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 + 1.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 1.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 1.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 2.0000000e+00 + 2.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 + 2.0000000e+00 2.0000000e+00 3.0000000e+00 3.0000000e+00 2.0000000e+00 1.0000000e+00 0.0000000e+00 2.0000000e+00 3.0000000e+00 0.0000000e+00 + 1.0000000e+00 3.0000000e+00 3.0000000e+00 0.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 + 0.0000000e+00 3.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 2.0000000e+00 3.0000000e+00 1.0000000e+00 1.0000000e+00 2.0000000e+00 + 2.0000000e+00 3.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 0.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 + 3.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 0.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 2.0000000e+00 0.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 3.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 + 0.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 1.0000000e+00 2.0000000e+00 2.0000000e+00 2.0000000e+00 1.0000000e+00 diff --git a/params/DbPlatform/initialTopographyRamp80m.txt b/params/DbPlatform/initialTopographyRamp80m.txt new file mode 100644 index 0000000..26d4887 --- /dev/null +++ b/params/DbPlatform/initialTopographyRamp80m.txt @@ -0,0 +1,150 @@ +1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 +1.53E+00 1.53E+00 1.53E+00 1.53E+00 1.53E+00 1.53E+00 1.53E+00 1.53E+00 1.53E+00 1.53E+00 +2.06E+00 2.06E+00 2.06E+00 2.06E+00 2.06E+00 2.06E+00 2.06E+00 2.06E+00 2.06E+00 2.06E+00 +2.59E+00 2.59E+00 2.59E+00 2.59E+00 2.59E+00 2.59E+00 2.59E+00 2.59E+00 2.59E+00 2.59E+00 +3.12E+00 3.12E+00 3.12E+00 3.12E+00 3.12E+00 3.12E+00 3.12E+00 3.12E+00 3.12E+00 3.12E+00 +3.65E+00 3.65E+00 3.65E+00 3.65E+00 3.65E+00 3.65E+00 3.65E+00 3.65E+00 3.65E+00 3.65E+00 +4.18E+00 4.18E+00 4.18E+00 4.18E+00 4.18E+00 4.18E+00 4.18E+00 4.18E+00 4.18E+00 4.18E+00 +4.71E+00 4.71E+00 4.71E+00 4.71E+00 4.71E+00 4.71E+00 4.71E+00 4.71E+00 4.71E+00 4.71E+00 +5.24E+00 5.24E+00 5.24E+00 5.24E+00 5.24E+00 5.24E+00 5.24E+00 5.24E+00 5.24E+00 5.24E+00 +5.77E+00 5.77E+00 5.77E+00 5.77E+00 5.77E+00 5.77E+00 5.77E+00 5.77E+00 5.77E+00 5.77E+00 +6.30E+00 6.30E+00 6.30E+00 6.30E+00 6.30E+00 6.30E+00 6.30E+00 6.30E+00 6.30E+00 6.30E+00 +6.83E+00 6.83E+00 6.83E+00 6.83E+00 6.83E+00 6.83E+00 6.83E+00 6.83E+00 6.83E+00 6.83E+00 +7.36E+00 7.36E+00 7.36E+00 7.36E+00 7.36E+00 7.36E+00 7.36E+00 7.36E+00 7.36E+00 7.36E+00 +7.89E+00 7.89E+00 7.89E+00 7.89E+00 7.89E+00 7.89E+00 7.89E+00 7.89E+00 7.89E+00 7.89E+00 +8.42E+00 8.42E+00 8.42E+00 8.42E+00 8.42E+00 8.42E+00 8.42E+00 8.42E+00 8.42E+00 8.42E+00 +8.95E+00 8.95E+00 8.95E+00 8.95E+00 8.95E+00 8.95E+00 8.95E+00 8.95E+00 8.95E+00 8.95E+00 +9.48E+00 9.48E+00 9.48E+00 9.48E+00 9.48E+00 9.48E+00 9.48E+00 9.48E+00 9.48E+00 9.48E+00 +1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+01 1.00E+01 +1.05E+01 1.05E+01 1.05E+01 1.05E+01 1.05E+01 1.05E+01 1.05E+01 1.05E+01 1.05E+01 1.05E+01 +1.11E+01 1.11E+01 1.11E+01 1.11E+01 1.11E+01 1.11E+01 1.11E+01 1.11E+01 1.11E+01 1.11E+01 +1.16E+01 1.16E+01 1.16E+01 1.16E+01 1.16E+01 1.16E+01 1.16E+01 1.16E+01 1.16E+01 1.16E+01 +1.21E+01 1.21E+01 1.21E+01 1.21E+01 1.21E+01 1.21E+01 1.21E+01 1.21E+01 1.21E+01 1.21E+01 +1.27E+01 1.27E+01 1.27E+01 1.27E+01 1.27E+01 1.27E+01 1.27E+01 1.27E+01 1.27E+01 1.27E+01 +1.32E+01 1.32E+01 1.32E+01 1.32E+01 1.32E+01 1.32E+01 1.32E+01 1.32E+01 1.32E+01 1.32E+01 +1.37E+01 1.37E+01 1.37E+01 1.37E+01 1.37E+01 1.37E+01 1.37E+01 1.37E+01 1.37E+01 1.37E+01 +1.43E+01 1.43E+01 1.43E+01 1.43E+01 1.43E+01 1.43E+01 1.43E+01 1.43E+01 1.43E+01 1.43E+01 +1.48E+01 1.48E+01 1.48E+01 1.48E+01 1.48E+01 1.48E+01 1.48E+01 1.48E+01 1.48E+01 1.48E+01 +1.53E+01 1.53E+01 1.53E+01 1.53E+01 1.53E+01 1.53E+01 1.53E+01 1.53E+01 1.53E+01 1.53E+01 +1.58E+01 1.58E+01 1.58E+01 1.58E+01 1.58E+01 1.58E+01 1.58E+01 1.58E+01 1.58E+01 1.58E+01 +1.64E+01 1.64E+01 1.64E+01 1.64E+01 1.64E+01 1.64E+01 1.64E+01 1.64E+01 1.64E+01 1.64E+01 +1.69E+01 1.69E+01 1.69E+01 1.69E+01 1.69E+01 1.69E+01 1.69E+01 1.69E+01 1.69E+01 1.69E+01 +1.74E+01 1.74E+01 1.74E+01 1.74E+01 1.74E+01 1.74E+01 1.74E+01 1.74E+01 1.74E+01 1.74E+01 +1.80E+01 1.80E+01 1.80E+01 1.80E+01 1.80E+01 1.80E+01 1.80E+01 1.80E+01 1.80E+01 1.80E+01 +1.85E+01 1.85E+01 1.85E+01 1.85E+01 1.85E+01 1.85E+01 1.85E+01 1.85E+01 1.85E+01 1.85E+01 +1.90E+01 1.90E+01 1.90E+01 1.90E+01 1.90E+01 1.90E+01 1.90E+01 1.90E+01 1.90E+01 1.90E+01 +1.96E+01 1.96E+01 1.96E+01 1.96E+01 1.96E+01 1.96E+01 1.96E+01 1.96E+01 1.96E+01 1.96E+01 +2.01E+01 2.01E+01 2.01E+01 2.01E+01 2.01E+01 2.01E+01 2.01E+01 2.01E+01 2.01E+01 2.01E+01 +2.06E+01 2.06E+01 2.06E+01 2.06E+01 2.06E+01 2.06E+01 2.06E+01 2.06E+01 2.06E+01 2.06E+01 +2.11E+01 2.11E+01 2.11E+01 2.11E+01 2.11E+01 2.11E+01 2.11E+01 2.11E+01 2.11E+01 2.11E+01 +2.17E+01 2.17E+01 2.17E+01 2.17E+01 2.17E+01 2.17E+01 2.17E+01 2.17E+01 2.17E+01 2.17E+01 +2.22E+01 2.22E+01 2.22E+01 2.22E+01 2.22E+01 2.22E+01 2.22E+01 2.22E+01 2.22E+01 2.22E+01 +2.27E+01 2.27E+01 2.27E+01 2.27E+01 2.27E+01 2.27E+01 2.27E+01 2.27E+01 2.27E+01 2.27E+01 +2.33E+01 2.33E+01 2.33E+01 2.33E+01 2.33E+01 2.33E+01 2.33E+01 2.33E+01 2.33E+01 2.33E+01 +2.38E+01 2.38E+01 2.38E+01 2.38E+01 2.38E+01 2.38E+01 2.38E+01 2.38E+01 2.38E+01 2.38E+01 +2.43E+01 2.43E+01 2.43E+01 2.43E+01 2.43E+01 2.43E+01 2.43E+01 2.43E+01 2.43E+01 2.43E+01 +2.49E+01 2.49E+01 2.49E+01 2.49E+01 2.49E+01 2.49E+01 2.49E+01 2.49E+01 2.49E+01 2.49E+01 +2.54E+01 2.54E+01 2.54E+01 2.54E+01 2.54E+01 2.54E+01 2.54E+01 2.54E+01 2.54E+01 2.54E+01 +2.59E+01 2.59E+01 2.59E+01 2.59E+01 2.59E+01 2.59E+01 2.59E+01 2.59E+01 2.59E+01 2.59E+01 +2.64E+01 2.64E+01 2.64E+01 2.64E+01 2.64E+01 2.64E+01 2.64E+01 2.64E+01 2.64E+01 2.64E+01 +2.70E+01 2.70E+01 2.70E+01 2.70E+01 2.70E+01 2.70E+01 2.70E+01 2.70E+01 2.70E+01 2.70E+01 +2.75E+01 2.75E+01 2.75E+01 2.75E+01 2.75E+01 2.75E+01 2.75E+01 2.75E+01 2.75E+01 2.75E+01 +2.80E+01 2.80E+01 2.80E+01 2.80E+01 2.80E+01 2.80E+01 2.80E+01 2.80E+01 2.80E+01 2.80E+01 +2.86E+01 2.86E+01 2.86E+01 2.86E+01 2.86E+01 2.86E+01 2.86E+01 2.86E+01 2.86E+01 2.86E+01 +2.91E+01 2.91E+01 2.91E+01 2.91E+01 2.91E+01 2.91E+01 2.91E+01 2.91E+01 2.91E+01 2.91E+01 +2.96E+01 2.96E+01 2.96E+01 2.96E+01 2.96E+01 2.96E+01 2.96E+01 2.96E+01 2.96E+01 2.96E+01 +3.02E+01 3.02E+01 3.02E+01 3.02E+01 3.02E+01 3.02E+01 3.02E+01 3.02E+01 3.02E+01 3.02E+01 +3.07E+01 3.07E+01 3.07E+01 3.07E+01 3.07E+01 3.07E+01 3.07E+01 3.07E+01 3.07E+01 3.07E+01 +3.12E+01 3.12E+01 3.12E+01 3.12E+01 3.12E+01 3.12E+01 3.12E+01 3.12E+01 3.12E+01 3.12E+01 +3.17E+01 3.17E+01 3.17E+01 3.17E+01 3.17E+01 3.17E+01 3.17E+01 3.17E+01 3.17E+01 3.17E+01 +3.23E+01 3.23E+01 3.23E+01 3.23E+01 3.23E+01 3.23E+01 3.23E+01 3.23E+01 3.23E+01 3.23E+01 +3.28E+01 3.28E+01 3.28E+01 3.28E+01 3.28E+01 3.28E+01 3.28E+01 3.28E+01 3.28E+01 3.28E+01 +3.33E+01 3.33E+01 3.33E+01 3.33E+01 3.33E+01 3.33E+01 3.33E+01 3.33E+01 3.33E+01 3.33E+01 +3.39E+01 3.39E+01 3.39E+01 3.39E+01 3.39E+01 3.39E+01 3.39E+01 3.39E+01 3.39E+01 3.39E+01 +3.44E+01 3.44E+01 3.44E+01 3.44E+01 3.44E+01 3.44E+01 3.44E+01 3.44E+01 3.44E+01 3.44E+01 +3.49E+01 3.49E+01 3.49E+01 3.49E+01 3.49E+01 3.49E+01 3.49E+01 3.49E+01 3.49E+01 3.49E+01 +3.55E+01 3.55E+01 3.55E+01 3.55E+01 3.55E+01 3.55E+01 3.55E+01 3.55E+01 3.55E+01 3.55E+01 +3.60E+01 3.60E+01 3.60E+01 3.60E+01 3.60E+01 3.60E+01 3.60E+01 3.60E+01 3.60E+01 3.60E+01 +3.65E+01 3.65E+01 3.65E+01 3.65E+01 3.65E+01 3.65E+01 3.65E+01 3.65E+01 3.65E+01 3.65E+01 +3.70E+01 3.70E+01 3.70E+01 3.70E+01 3.70E+01 3.70E+01 3.70E+01 3.70E+01 3.70E+01 3.70E+01 +3.76E+01 3.76E+01 3.76E+01 3.76E+01 3.76E+01 3.76E+01 3.76E+01 3.76E+01 3.76E+01 3.76E+01 +3.81E+01 3.81E+01 3.81E+01 3.81E+01 3.81E+01 3.81E+01 3.81E+01 3.81E+01 3.81E+01 3.81E+01 +3.86E+01 3.86E+01 3.86E+01 3.86E+01 3.86E+01 3.86E+01 3.86E+01 3.86E+01 3.86E+01 3.86E+01 +3.92E+01 3.92E+01 3.92E+01 3.92E+01 3.92E+01 3.92E+01 3.92E+01 3.92E+01 3.92E+01 3.92E+01 +3.97E+01 3.97E+01 3.97E+01 3.97E+01 3.97E+01 3.97E+01 3.97E+01 3.97E+01 3.97E+01 3.97E+01 +4.02E+01 4.02E+01 4.02E+01 4.02E+01 4.02E+01 4.02E+01 4.02E+01 4.02E+01 4.02E+01 4.02E+01 +4.08E+01 4.08E+01 4.08E+01 4.08E+01 4.08E+01 4.08E+01 4.08E+01 4.08E+01 4.08E+01 4.08E+01 +4.13E+01 4.13E+01 4.13E+01 4.13E+01 4.13E+01 4.13E+01 4.13E+01 4.13E+01 4.13E+01 4.13E+01 +4.18E+01 4.18E+01 4.18E+01 4.18E+01 4.18E+01 4.18E+01 4.18E+01 4.18E+01 4.18E+01 4.18E+01 +4.23E+01 4.23E+01 4.23E+01 4.23E+01 4.23E+01 4.23E+01 4.23E+01 4.23E+01 4.23E+01 4.23E+01 +4.29E+01 4.29E+01 4.29E+01 4.29E+01 4.29E+01 4.29E+01 4.29E+01 4.29E+01 4.29E+01 4.29E+01 +4.34E+01 4.34E+01 4.34E+01 4.34E+01 4.34E+01 4.34E+01 4.34E+01 4.34E+01 4.34E+01 4.34E+01 +4.39E+01 4.39E+01 4.39E+01 4.39E+01 4.39E+01 4.39E+01 4.39E+01 4.39E+01 4.39E+01 4.39E+01 +4.45E+01 4.45E+01 4.45E+01 4.45E+01 4.45E+01 4.45E+01 4.45E+01 4.45E+01 4.45E+01 4.45E+01 +4.50E+01 4.50E+01 4.50E+01 4.50E+01 4.50E+01 4.50E+01 4.50E+01 4.50E+01 4.50E+01 4.50E+01 +4.55E+01 4.55E+01 4.55E+01 4.55E+01 4.55E+01 4.55E+01 4.55E+01 4.55E+01 4.55E+01 4.55E+01 +4.61E+01 4.61E+01 4.61E+01 4.61E+01 4.61E+01 4.61E+01 4.61E+01 4.61E+01 4.61E+01 4.61E+01 +4.66E+01 4.66E+01 4.66E+01 4.66E+01 4.66E+01 4.66E+01 4.66E+01 4.66E+01 4.66E+01 4.66E+01 +4.71E+01 4.71E+01 4.71E+01 4.71E+01 4.71E+01 4.71E+01 4.71E+01 4.71E+01 4.71E+01 4.71E+01 +4.76E+01 4.76E+01 4.76E+01 4.76E+01 4.76E+01 4.76E+01 4.76E+01 4.76E+01 4.76E+01 4.76E+01 +4.82E+01 4.82E+01 4.82E+01 4.82E+01 4.82E+01 4.82E+01 4.82E+01 4.82E+01 4.82E+01 4.82E+01 +4.87E+01 4.87E+01 4.87E+01 4.87E+01 4.87E+01 4.87E+01 4.87E+01 4.87E+01 4.87E+01 4.87E+01 +4.92E+01 4.92E+01 4.92E+01 4.92E+01 4.92E+01 4.92E+01 4.92E+01 4.92E+01 4.92E+01 4.92E+01 +4.98E+01 4.98E+01 4.98E+01 4.98E+01 4.98E+01 4.98E+01 4.98E+01 4.98E+01 4.98E+01 4.98E+01 +5.03E+01 5.03E+01 5.03E+01 5.03E+01 5.03E+01 5.03E+01 5.03E+01 5.03E+01 5.03E+01 5.03E+01 +5.08E+01 5.08E+01 5.08E+01 5.08E+01 5.08E+01 5.08E+01 5.08E+01 5.08E+01 5.08E+01 5.08E+01 +5.14E+01 5.14E+01 5.14E+01 5.14E+01 5.14E+01 5.14E+01 5.14E+01 5.14E+01 5.14E+01 5.14E+01 +5.19E+01 5.19E+01 5.19E+01 5.19E+01 5.19E+01 5.19E+01 5.19E+01 5.19E+01 5.19E+01 5.19E+01 +5.24E+01 5.24E+01 5.24E+01 5.24E+01 5.24E+01 5.24E+01 5.24E+01 5.24E+01 5.24E+01 5.24E+01 +5.29E+01 5.29E+01 5.29E+01 5.29E+01 5.29E+01 5.29E+01 5.29E+01 5.29E+01 5.29E+01 5.29E+01 +5.35E+01 5.35E+01 5.35E+01 5.35E+01 5.35E+01 5.35E+01 5.35E+01 5.35E+01 5.35E+01 5.35E+01 +5.40E+01 5.40E+01 5.40E+01 5.40E+01 5.40E+01 5.40E+01 5.40E+01 5.40E+01 5.40E+01 5.40E+01 +5.45E+01 5.45E+01 5.45E+01 5.45E+01 5.45E+01 5.45E+01 5.45E+01 5.45E+01 5.45E+01 5.45E+01 +5.51E+01 5.51E+01 5.51E+01 5.51E+01 5.51E+01 5.51E+01 5.51E+01 5.51E+01 5.51E+01 5.51E+01 +5.56E+01 5.56E+01 5.56E+01 5.56E+01 5.56E+01 5.56E+01 5.56E+01 5.56E+01 5.56E+01 5.56E+01 +5.61E+01 5.61E+01 5.61E+01 5.61E+01 5.61E+01 5.61E+01 5.61E+01 5.61E+01 5.61E+01 5.61E+01 +5.67E+01 5.67E+01 5.67E+01 5.67E+01 5.67E+01 5.67E+01 5.67E+01 5.67E+01 5.67E+01 5.67E+01 +5.72E+01 5.72E+01 5.72E+01 5.72E+01 5.72E+01 5.72E+01 5.72E+01 5.72E+01 5.72E+01 5.72E+01 +5.77E+01 5.77E+01 5.77E+01 5.77E+01 5.77E+01 5.77E+01 5.77E+01 5.77E+01 5.77E+01 5.77E+01 +5.82E+01 5.82E+01 5.82E+01 5.82E+01 5.82E+01 5.82E+01 5.82E+01 5.82E+01 5.82E+01 5.82E+01 +5.88E+01 5.88E+01 5.88E+01 5.88E+01 5.88E+01 5.88E+01 5.88E+01 5.88E+01 5.88E+01 5.88E+01 +5.93E+01 5.93E+01 5.93E+01 5.93E+01 5.93E+01 5.93E+01 5.93E+01 5.93E+01 5.93E+01 5.93E+01 +5.98E+01 5.98E+01 5.98E+01 5.98E+01 5.98E+01 5.98E+01 5.98E+01 5.98E+01 5.98E+01 5.98E+01 +6.04E+01 6.04E+01 6.04E+01 6.04E+01 6.04E+01 6.04E+01 6.04E+01 6.04E+01 6.04E+01 6.04E+01 +6.09E+01 6.09E+01 6.09E+01 6.09E+01 6.09E+01 6.09E+01 6.09E+01 6.09E+01 6.09E+01 6.09E+01 +6.14E+01 6.14E+01 6.14E+01 6.14E+01 6.14E+01 6.14E+01 6.14E+01 6.14E+01 6.14E+01 6.14E+01 +6.20E+01 6.20E+01 6.20E+01 6.20E+01 6.20E+01 6.20E+01 6.20E+01 6.20E+01 6.20E+01 6.20E+01 +6.25E+01 6.25E+01 6.25E+01 6.25E+01 6.25E+01 6.25E+01 6.25E+01 6.25E+01 6.25E+01 6.25E+01 +6.30E+01 6.30E+01 6.30E+01 6.30E+01 6.30E+01 6.30E+01 6.30E+01 6.30E+01 6.30E+01 6.30E+01 +6.35E+01 6.35E+01 6.35E+01 6.35E+01 6.35E+01 6.35E+01 6.35E+01 6.35E+01 6.35E+01 6.35E+01 +6.41E+01 6.41E+01 6.41E+01 6.41E+01 6.41E+01 6.41E+01 6.41E+01 6.41E+01 6.41E+01 6.41E+01 +6.46E+01 6.46E+01 6.46E+01 6.46E+01 6.46E+01 6.46E+01 6.46E+01 6.46E+01 6.46E+01 6.46E+01 +6.51E+01 6.51E+01 6.51E+01 6.51E+01 6.51E+01 6.51E+01 6.51E+01 6.51E+01 6.51E+01 6.51E+01 +6.57E+01 6.57E+01 6.57E+01 6.57E+01 6.57E+01 6.57E+01 6.57E+01 6.57E+01 6.57E+01 6.57E+01 +6.62E+01 6.62E+01 6.62E+01 6.62E+01 6.62E+01 6.62E+01 6.62E+01 6.62E+01 6.62E+01 6.62E+01 +6.67E+01 6.67E+01 6.67E+01 6.67E+01 6.67E+01 6.67E+01 6.67E+01 6.67E+01 6.67E+01 6.67E+01 +6.73E+01 6.73E+01 6.73E+01 6.73E+01 6.73E+01 6.73E+01 6.73E+01 6.73E+01 6.73E+01 6.73E+01 +6.78E+01 6.78E+01 6.78E+01 6.78E+01 6.78E+01 6.78E+01 6.78E+01 6.78E+01 6.78E+01 6.78E+01 +6.83E+01 6.83E+01 6.83E+01 6.83E+01 6.83E+01 6.83E+01 6.83E+01 6.83E+01 6.83E+01 6.83E+01 +6.88E+01 6.88E+01 6.88E+01 6.88E+01 6.88E+01 6.88E+01 6.88E+01 6.88E+01 6.88E+01 6.88E+01 +6.94E+01 6.94E+01 6.94E+01 6.94E+01 6.94E+01 6.94E+01 6.94E+01 6.94E+01 6.94E+01 6.94E+01 +6.99E+01 6.99E+01 6.99E+01 6.99E+01 6.99E+01 6.99E+01 6.99E+01 6.99E+01 6.99E+01 6.99E+01 +7.04E+01 7.04E+01 7.04E+01 7.04E+01 7.04E+01 7.04E+01 7.04E+01 7.04E+01 7.04E+01 7.04E+01 +7.10E+01 7.10E+01 7.10E+01 7.10E+01 7.10E+01 7.10E+01 7.10E+01 7.10E+01 7.10E+01 7.10E+01 +7.15E+01 7.15E+01 7.15E+01 7.15E+01 7.15E+01 7.15E+01 7.15E+01 7.15E+01 7.15E+01 7.15E+01 +7.20E+01 7.20E+01 7.20E+01 7.20E+01 7.20E+01 7.20E+01 7.20E+01 7.20E+01 7.20E+01 7.20E+01 +7.26E+01 7.26E+01 7.26E+01 7.26E+01 7.26E+01 7.26E+01 7.26E+01 7.26E+01 7.26E+01 7.26E+01 +7.31E+01 7.31E+01 7.31E+01 7.31E+01 7.31E+01 7.31E+01 7.31E+01 7.31E+01 7.31E+01 7.31E+01 +7.36E+01 7.36E+01 7.36E+01 7.36E+01 7.36E+01 7.36E+01 7.36E+01 7.36E+01 7.36E+01 7.36E+01 +7.41E+01 7.41E+01 7.41E+01 7.41E+01 7.41E+01 7.41E+01 7.41E+01 7.41E+01 7.41E+01 7.41E+01 +7.47E+01 7.47E+01 7.47E+01 7.47E+01 7.47E+01 7.47E+01 7.47E+01 7.47E+01 7.47E+01 7.47E+01 +7.52E+01 7.52E+01 7.52E+01 7.52E+01 7.52E+01 7.52E+01 7.52E+01 7.52E+01 7.52E+01 7.52E+01 +7.57E+01 7.57E+01 7.57E+01 7.57E+01 7.57E+01 7.57E+01 7.57E+01 7.57E+01 7.57E+01 7.57E+01 +7.63E+01 7.63E+01 7.63E+01 7.63E+01 7.63E+01 7.63E+01 7.63E+01 7.63E+01 7.63E+01 7.63E+01 +7.68E+01 7.68E+01 7.68E+01 7.68E+01 7.68E+01 7.68E+01 7.68E+01 7.68E+01 7.68E+01 7.68E+01 +7.73E+01 7.73E+01 7.73E+01 7.73E+01 7.73E+01 7.73E+01 7.73E+01 7.73E+01 7.73E+01 7.73E+01 +7.79E+01 7.79E+01 7.79E+01 7.79E+01 7.79E+01 7.79E+01 7.79E+01 7.79E+01 7.79E+01 7.79E+01 +7.84E+01 7.84E+01 7.84E+01 7.84E+01 7.84E+01 7.84E+01 7.84E+01 7.84E+01 7.84E+01 7.84E+01 +7.89E+01 7.89E+01 7.89E+01 7.89E+01 7.89E+01 7.89E+01 7.89E+01 7.89E+01 7.89E+01 7.89E+01 +7.94E+01 7.94E+01 7.94E+01 7.94E+01 7.94E+01 7.94E+01 7.94E+01 7.94E+01 7.94E+01 7.94E+01 +8.00E+01 8.00E+01 8.00E+01 8.00E+01 8.00E+01 8.00E+01 8.00E+01 8.00E+01 8.00E+01 8.00E+01 \ No newline at end of file diff --git a/params/DbPlatform/paramsConcentration.txt b/params/DbPlatform/paramsConcentration.txt new file mode 100644 index 0000000..830debf --- /dev/null +++ b/params/DbPlatform/paramsConcentration.txt @@ -0,0 +1,8 @@ +10000000 carbonate input at source point(in m3/yrs) +10 y position of input carbonate +25 x position of input carbonate +1000000 diffusion on y direction plus, in m2/yr +500000 diffusion on y direction minus, in m2/yr +500000 diffusion on x direction plus, in m2/yr +100000 diffusion on x direction minus, in m2/yr +0.0 initial concentration value for the whole grid, in kg/m3 diff --git a/params/DbPlatform/paramsInputValues.txt b/params/DbPlatform/paramsInputValues.txt new file mode 100644 index 0000000..db3a216 --- /dev/null +++ b/params/DbPlatform/paramsInputValues.txt @@ -0,0 +1,64 @@ +LongRamp Model Name +10 Total number of grid cells in the x direction (Strike) +150 Total number of grid cells in the y direction (dip) +100 Grid cell edge dimension, in meters +200 Total iterations +0.01 Timestep( My) + +1.00 Eustasy period 1 (My) +20.0 Eustasy amplitude 1 (m) +0.112 Eustasy period 2 (My) +5.0 Eustasy amplitude 2 (m) + +1 Timesteps per CA iteration at max prod rate +1 Timesteps per CA iteration at min prod rate + +100000 Bathymetry difference between cells that get disconected in CA + +3 Number of producing factories + +500 Lithology 1 REEF maximum carb production Rate (mMy-1) +0.3 Lithology 1 Surface light intensity / bell shape, depth of maximum production +30 Lithology 1 extinction coefficient / bell shape, vertical width of the bell +20 Lithology 1 saturating light / bell shape, shape of the curve, rounded or spiky +4 Lithology 1 transport product facies +0.9 Lithology 1 transported fraction of total production +0.3 Lithology 1 minimum gradient for transport, m per m +0.7 Lithology 1 probability threshold for gradient-independent continuous transport +0.9 Minimum tolerated wave energy, reef +1.0 Maximum tolerated wave energy, reef + +160 Lithology 2 Platform interior maximum carb production Rate (mMy-1) +0.1 Lithology 2 Surface light intensity / bell shape, depth of maximum production +40 Lithology 2 extinction coefficient / bell shape, vertical width of the bell +5 Lithology 2 saturating light / bell shape, shape of the curve, rounded or spiky +5 Lithology 2 transport product facies +0.4 Lithology 2 transported fraction of total production +0.2 Lithology 2 minimum gradient for transport, m per m +0.7 Lithology 2 probability threshold for gradient-independent continuous transport +0.2 Minimum tolerated wave energy, platform interior +0.9 Maximum tolerated wave energy, platform interior + +150 Lithology 3 Platform interior maximum carb production Rate (mMy-1) +0.1 Lithology 3 Surface light intensity / bell shape, depth of maximum production +40 Lithology 3 extinction coefficient / bell shape, vertical width of the bell +10 Lithology 3 saturating light / bell shape, shape of the curve, rounded or spiky +6 Lithology 3 transport product facies +0.4 Lithology 3 transported fraction of total production +0.1 Lithology 3 minimum gradient for transport, m per m +0.7 Lithology 3 probability threshold for gradient-independent continuous transport +0.0 Minimum tolerated wave energy, platform interior +0.2 Maximum tolerated wave energy, platform interior + +params/DbPlatform/simpleOneRateSubsidence70mY150X10.txt spatially uniform constant subsdience rate m/My +params/DbPlatform/CARules3_4_10_6_10.txt CA rules input file +params/DbPlatform/initialFaciesMapRandomY150X10.txt initial condition facies map +params/DbPlatform/initialTopographyRamp80m.txt initial water depth map +params/DbPlatform/paramsConcentration.txt concentration input parameters +params/DbPlatform/paramsSiliciclastic.txt silicilastic input parameters +params/DbPlatform/seaLevelConstS1000iterations.txt sea level curve, maybe not needed if sinusoid selected +params/DbPlatform/carbProdConst1000Iterations.txt carbonate maximum production rate time curve +colorMaps/carboCAT3FaciesColours.txt facies colour map files + +10 Number of time lines to plot on cross sections +50 100 150 200 250 300 350 400 450 500 Ages of the time lines in terms of model iteration number diff --git a/params/DbPlatform/paramsProcesses.txt b/params/DbPlatform/paramsProcesses.txt new file mode 100644 index 0000000..40fb533 --- /dev/null +++ b/params/DbPlatform/paramsProcesses.txt @@ -0,0 +1,9 @@ +simpleWave % Facies distribution routines: 'simpleWave','orderedCA' or 'productionCA' or 'waveCA' or 'waveHybridCA' +steepestSlope % transportation algorithm: 'steepestSlope' or 'crossPlatform' +on % wave Energy +off % siliciciclastics: 'on' or 'off' +off % carbonate concentration: 'on' or 'off' +off % soil deposition: 'on' or 'off' +pre-exposed % cell state when re-flooding: 'empty' or 'pre-exposed' +sinusoid % sea level curve origin: 'sinusoid' or 'file' +wrap % how to check neighbours or transport: 'wrap' or 'unwrap' \ No newline at end of file diff --git a/params/DbPlatform/paramsSiliciclastic.txt b/params/DbPlatform/paramsSiliciclastic.txt new file mode 100644 index 0000000..fb5e279 --- /dev/null +++ b/params/DbPlatform/paramsSiliciclastic.txt @@ -0,0 +1,6 @@ +50 siliciclastic at source point(in m) +1 y position of input carbonate +25 x position of input carbonate +0 source length (0 for point source, >0 for line source) +2000 diffusion on y direction, in m2/yr +2000 diffusion on x direction, in m2/yr \ No newline at end of file diff --git a/params/DbPlatform/seaLevelConst1000iterations.txt b/params/DbPlatform/seaLevelConst1000iterations.txt new file mode 100644 index 0000000..d7d54a3 --- /dev/null +++ b/params/DbPlatform/seaLevelConst1000iterations.txt @@ -0,0 +1,1001 @@ +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/params/DbPlatform/simpleOneRateSubsidence70mY150X10.txt b/params/DbPlatform/simpleOneRateSubsidence70mY150X10.txt new file mode 100644 index 0000000..e55447b --- /dev/null +++ b/params/DbPlatform/simpleOneRateSubsidence70mY150X10.txt @@ -0,0 +1,150 @@ +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 +7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 7.00E+01 \ No newline at end of file diff --git a/plot3DBlock.m b/plot3DBlock.m new file mode 100644 index 0000000..0fd4c00 --- /dev/null +++ b/plot3DBlock.m @@ -0,0 +1,75 @@ +function [glob,graph] = plot3DBlock(glob, iteration) + + % ScreenSize is a four-element vector: [left, bottom, width, height]: + scrsz = get(0,'ScreenSize'); % vector + % Position requires left bottom width height values. screensize vector is in this format 1=left 2=bottom 3=width 4=height + figure('Visible','on','Position',[1 scrsz(2)+20 scrsz(3)*0.8 scrsz(4)*0.8]); + + plotCrossSectionDip3D(1, glob, iteration); + plotCrossSectionStrike3D(glob.ySize, glob, iteration); + plot3DPlanView(glob, iteration); + + % Draw the final sea-level as a rectangle + yco = [glob.dx, glob.ySize*glob.dx, glob.ySize*glob.dx, glob.dx]; + zco = [glob.SL(iteration), glob.SL(iteration), glob.SL(iteration), glob.SL(iteration)]; + xco = [glob.dx, glob.dx, (glob.xSize+1)*glob.dx, (glob.xSize+1)*glob.dx]; + line(xco, yco, zco, 'LineWidth',2, 'color', 'blue'); + + xlabel('X Distance (m)','FontSize',11); + ylabel('Y Distance (m)','FontSize',11); + zlabel('Elevation (m)','FontSize',11); + grid on; + view(-135,50); + ax = gca; + ax.DataAspectRatio = [1, 3, 0.03]; % x-scale=1, y-scale = 3, z-scale = 0.002 Change these scalings to alter the aspect ratio of the 3d plot +end + +function plot3DPlanView(glob, iteration) + + % Loop along the section and through timelines to draw the cross section + for x = 1:glob.xSize + + for y = 1:glob.ySize-1 + + % Find the youngest iteration that has deposition + t = iteration; + while ((glob.numberOfLayers(y,x, t)) == 0 || sum(glob.thickness{y,x,t}(:)) < 0.01) && t > 1 + t = t - 1; + end + + topCell = glob.numberOfLayers(y,x, t); + if topCell > 0 + + fCode = glob.facies{y,x,t}(topCell); % Note this is zero for no depositon, 9 for subaerial hiatus + if fCode > 0 + faciesCol = [glob.faciesColours(fCode,1) glob.faciesColours(fCode,2) glob.faciesColours(fCode,3)]; + else + faciesCol = [1,1,1]; + end + + % Draw the top planform patch + zco = [glob.strata(y,x,t), glob.strata(y,x,t), glob.strata(y,x,t), glob.strata(y,x,t)]; + xco = [((x-0.5)*glob.dx), ((x+0.5)*glob.dx), ((x+0.5)*glob.dx), ((x-0.5)*glob.dx)]; + yco = [double(y-0.5)*glob.dx, double(y-0.5)*glob.dx, double(y+0.5)*glob.dx, double(y+0.5)*glob.dx]; + patch(xco, yco, zco, faciesCol , 'EdgeColor','none'); + + % Draw the front vertical panel patch + zco = [glob.strata(y,x,t), glob.strata(y,x,t), glob.strata(y+1,x,t), glob.strata(y+1,x,t)]; + xco = [((x-0.5)*glob.dx), ((x+0.5)*glob.dx), ((x+0.5)*glob.dx), ((x-0.5)*glob.dx)]; + yco = [double(y+0.5)*glob.dx, double(y+0.5)*glob.dx, double(y+0.5)*glob.dx, double(y+0.5)*glob.dx]; + patch(xco, yco, zco, faciesCol , 'EdgeColor','none'); + end + end + end + + % Draw dip-direction line along cell boundaries + for x = 1:glob.xSize + xco = ones(1,glob.ySize) .* double(x*glob.dx); + yco = double(1:glob.ySize) .* glob.dx; + zco = glob.strata(1:glob.ySize,x,iteration); + line(xco, yco, zco, 'LineWidth',2, 'color', 'black'); + end + + xco = ones(1,glob.ySize) .* double((x+1)*glob.dx); + line(xco, yco, zco, 'LineWidth',2, 'color', 'black'); +end \ No newline at end of file diff --git a/plotChronostratSectionDip.m b/plotChronostratSectionDip.m new file mode 100644 index 0000000..6bca9f7 --- /dev/null +++ b/plotChronostratSectionDip.m @@ -0,0 +1,47 @@ +function plotChronostratSectionDip( xPos,chronoPlot, glob, iteration) + + subplot(chronoPlot); + thicknessPerY=sum(glob.numberOfLayers(:,xPos,2:end),3); + longA=find(thicknessPerY>0, 1, 'last' ); + + for i=0:10 + patch(0,0,i); % use a dummy patch to force colour map to range 0-10 + end + + for k=1:iteration-1 + for y = 1:longA + % xco is a vector containing x-axis coords for the corner points of + % each strat section grid cell. Yco is the equivalent y-axis vector + xco = [y*glob.dx, y*glob.dx, (y+1)*glob.dx,(y+1)*glob.dx]; + + faciesList = glob.facies{y,xPos,k+1}; + oneThickness = glob.strata(y,xPos, k+1) - glob.strata(y,xPos,k); + + if max(faciesList) > 0 + cellHeight = glob.deltaT / (glob.numberOfLayers(y,xPos,k+1)); + else + cellHeight = glob.deltaT; + end + cellBase = k*glob.deltaT; + + % Now draw the facies, however many there are... + % cellBase = cellBase+cellHeight; % to account for in-situ facies cell + for fLoop = 1:glob.numberOfLayers(y,xPos,k+1) + + tco = [cellBase, cellBase+cellHeight, cellBase+cellHeight, cellBase]; + fCode = faciesList(fLoop); + if fCode > 0 && oneThickness > 0.01 + faciesCol = [glob.faciesColours(fCode,1) glob.faciesColours(fCode,2) glob.faciesColours(fCode,3)]; + patch(xco, tco, faciesCol,'EdgeColor','none'); + end + cellBase = cellBase + cellHeight; + end + end + end + + % Force 2 ticks on the y axis + axis([1*glob.dx, (longA+1)*glob.dx, 0, glob.deltaT * iteration]); + xlabel('Distance (km)','FontSize',11); + ylabel('E.M.T. (My)','FontSize',11); + set(gca,'FontSize',11) +end \ No newline at end of file diff --git a/plotCrossSectionDip2D.m b/plotCrossSectionDip2D.m new file mode 100644 index 0000000..715ba8c --- /dev/null +++ b/plotCrossSectionDip2D.m @@ -0,0 +1,82 @@ +function plotCrossSectionDip2D(x, crossSectionPlot, glob, iteration) + + subplot(crossSectionPlot); + thicknessPerY=sum(glob.numberOfLayers(:,x,2:end),3); + long = find(thicknessPerY>0, 1, 'last' ); + if isempty(long) + long=glob.ySize; + end + + % First delete the previous plot with a single large white patch + minDepth = max(max(glob.strata(:,:,iteration))); % Find the highest (so shallowest) values in the strata array + minDepth = minDepth * 1.1; % Add 10% to the minimum depth + maxDepth = min(min(glob.strata(:,:,1))); % Find the lowest (ie deepest) values in strata array + maxDepth = maxDepth * 1.1; % Add 10% to the maximum depth + + patch([0 long*glob.dx long*glob.dx 0], [minDepth minDepth maxDepth maxDepth], [1 1 1], 'EdgeColor','none'); + + % Now reuse maxdepth to draw a light grey solid colour basement at the bottom of the section + for y = 1:long%glob.ySize-1 + xco = [y*glob.dx, y*glob.dx, (y+1)*glob.dx, (y+1)*glob.dx]; + zco = [maxDepth, glob.strata(y,x,1), glob.strata(y,x,1), maxDepth]; + patch(xco, zco, [0.7 0.7 0.7],'EdgeColor','none'); + end + + % Loop along the section and through timelines to draw the cross section + for y = 1:long + for k=2:iteration + + cell = glob.numberOfLayers(y,x,k); + while cell > 0 + + xco = [y*glob.dx, y*glob.dx, (y+1)*glob.dx, (y+1)*glob.dx]; + + allThick = sum(glob.thickness{y,x,k}(1:cell)); + oneThick = sum(glob.thickness{y,x,k}(1:cell-1)); + top=glob.strata(y,x,k-1)+allThick; + bottom=glob.strata(y,x,k-1)+oneThick; + fCode = glob.facies{y,x,k}(cell); % Note this is zero for no depositon, 9 for subaerial hiatus + + % Draw the insitu production facies first + if fCode > 0 + zco = [bottom, top,top, bottom]; + faciesCol=[glob.faciesColours(fCode,1) glob.faciesColours(fCode,2) glob.faciesColours(fCode,3)] ; + patch(xco, zco, faciesCol,'EdgeColor','none'); + end + + cell = cell - 1; + end + end + end + + % Loop through iterations and draw timelines + for i=1:glob.timeLineCount + + k = glob.timeLineAge(i); + + if k <= iteration + for y = 1:long-1%glob.ySize-1 + % draw a marker line across the top and down/up the side of a particular grid cell + xco = [y*glob.dx, (y+1)*glob.dx, (y+1)*glob.dx]; + yco = [glob.strata(y,x,k), glob.strata(y,x,k), glob.strata(y+1,x,k)]; + line(xco, yco, 'LineWidth',2, 'color', 'black'); + end + line([(y+1)*glob.dx,(y+2)*glob.dx],[glob.strata(y+1,x,k), glob.strata(y+1,x,k)],'LineWidth',2, 'color', 'black'); + end + end + + % Draw the final sea-level + xco = [1*glob.dx, (long+1)*glob.dx]; + yco = [glob.SL(iteration) glob.SL(iteration)]; + line(xco,yco, 'LineWidth',2, 'color', 'blue'); + + if minDepth>glob.SL(iteration); mm=minDepth; else mm=glob.SL(iteration); end + axis([1 (long+1)*glob.dx maxDepth mm]); + l=numel(num2str(round(maxDepth)))-1; + maxD=(10^(l-1))*(ceil(maxDepth/(10^(l-1)))); + + ylabel('Elevation (m)','FontSize',11); + set(gca,'FontSize',11) + set(crossSectionPlot,'XTickLabel',[]); + grid on; +end \ No newline at end of file diff --git a/plotCrossSectionDip3D.m b/plotCrossSectionDip3D.m new file mode 100644 index 0000000..0ef3a9a --- /dev/null +++ b/plotCrossSectionDip3D.m @@ -0,0 +1,76 @@ +function plotCrossSectionDip3D(x, glob, iteration) + + thicknessPerY=sum(glob.numberOfLayers(:,x,2:end),3); + long = find(thicknessPerY>0, 1, 'last' ); + if isempty(long) + long=glob.ySize; + end + + % First find the plot limits for the z-axis + minDepth = max(max(glob.strata(:,:,iteration))); % Find the highest (so shallowest) values in the strata array + minDepth = minDepth * 1.1; % Add 10% to the minimum depth + maxDepth = min(min(glob.strata(:,:,1))); % Find the lowest (ie deepest) values in strata array + maxDepth = maxDepth * 1.1; % Add 10% to the maximum depth + +% % Use maxdepth to draw a light-grey solid colour basement along the x-axis bottom of the section +% % Use j as the loop counter, because x is used as the line of section x-coord +% yco = [double(glob.ySize)*glob.dx, double(glob.ySize)*glob.dx, double(glob.ySize)*glob.dx, double(glob.ySize)*glob.dx]; +% for j = 1:glob.xSize +% xco = [double(j-0.5)*glob.dx, double(j-0.5)*glob.dx, double(j+0.5)*glob.dx, double(j+0.5)*glob.dx]; +% zco = [maxDepth, glob.strata(glob.ySize,j,1), glob.strata(glob.ySize,j,1), maxDepth]; +% patch(xco, yco, zco, [0.7 0.7 0.7],'EdgeColor','none'); +% end + + % Because everything else will plot along the cross-section line at x, set xco here for use below + xco = [x*glob.dx,x*glob.dx,x*glob.dx,x*glob.dx]; % Everything will plot along cross-section line x=1 + + % Now use maxdepth to draw a light grey solid colour basement along the y-axis bottom of the section + for y = 1:long%glob.ySize-1 + yco = [double(y-0.5)*glob.dx, double(y-0.5)*glob.dx, double(y+0.5)*glob.dx, double(y+0.5)*glob.dx]; + zco = [maxDepth, glob.strata(y,x,1), glob.strata(y,x,1), maxDepth]; + patch(xco, yco, zco, [0.7 0.7 0.7],'EdgeColor','none'); + end + + % Loop along the section and through timelines to draw the cross section + for y = 1:long + for k=2:iteration + + cell = glob.numberOfLayers(y,x,k); + while cell > 0 + + yco = [double(y-0.5)*glob.dx, double(y-0.5)*glob.dx, double(y+0.5)*glob.dx, double(y+0.5)*glob.dx]; + + allThick = sum(glob.thickness{y,x,k}(1:cell)); + oneThick = sum(glob.thickness{y,x,k}(1:cell-1)); + top = glob.strata(y,x,k-1)+allThick; + bottom = glob.strata(y,x,k-1)+oneThick; + fCode = glob.facies{y,x,k}(cell); % Note this is zero for no depositon, 9 for subaerial hiatus + + % Draw the in-situ production facies first + if fCode > 0 + zco = [bottom, top, top, bottom]; + faciesCol = [glob.faciesColours(fCode,1) glob.faciesColours(fCode,2) glob.faciesColours(fCode,3)]; + patch(xco, yco, zco, faciesCol,'EdgeColor','none'); + end + + cell = cell - 1; + end + end + end + + % Loop through iterations and draw timelines + for i=1:glob.timeLineCount + + k = glob.timeLineAge(i); + + if k <= iteration + for y = 1:long-1%glob.ySize-1 + % Draw a marker line across the top and down/up the side of a particular grid cell + yco = [y*glob.dx, (y+1)*glob.dx, (y+1)*glob.dx]; + zco = [glob.strata(y,x,k), glob.strata(y,x,k), glob.strata(y+1,x,k)]; + line([x*glob.dx,x*glob.dx,x*glob.dx], yco, zco, 'LineWidth',2, 'color', 'black'); + end + line([x*glob.dx,x*glob.dx], [(y+1)*glob.dx,(y+2)*glob.dx], [glob.strata(y+1,x,k), glob.strata(y+1,x,k)],'LineWidth',2, 'color', 'black'); + end + end +end \ No newline at end of file diff --git a/plotCrossSectionStrike3D.m b/plotCrossSectionStrike3D.m new file mode 100644 index 0000000..0bf0ead --- /dev/null +++ b/plotCrossSectionStrike3D.m @@ -0,0 +1,67 @@ +function plotCrossSectionStrike3D(y, glob, iteration) + + thicknessPerX=sum(glob.numberOfLayers(y,:,2:end),3); + long = find(thicknessPerX>0, 1, 'last' ); + if isempty(long) + long=glob.xSize; + end + + % First find the plot limits for the z-axis + minDepth = max(max(glob.strata(:,:,iteration))); % Find the highest (so shallowest) values in the strata array + minDepth = minDepth * 1.1; % Add 10% to the minimum depth + maxDepth = min(min(glob.strata(:,:,1))); % Find the lowest (ie deepest) values in strata array + maxDepth = maxDepth * 1.1; % Add 10% to the maximum depth + + % Because everything plots along the cross-section line at y, set yco here for use below + yco = [double(y)*glob.dx,double(y)*glob.dx,double(y)*glob.dx,double(y)*glob.dx]; % Everything will plot along cross-section line x=1 + + % Now use maxdepth to draw a light grey solid colour basement along the y-axis bottom of the section + for x = 1:long%glob.ySize-1 + xco = [double(x-0.5)*glob.dx, double(x-0.5)*glob.dx, double(x+0.5)*glob.dx, double(x+0.5)*glob.dx]; + zco = [maxDepth, glob.strata(y,x,1), glob.strata(y,x,1), maxDepth]; + patch(xco, yco, zco, [0.7 0.7 0.7],'EdgeColor','none'); + end + + % Loop along the section and through timelines to draw the cross section + for x = 1:long + for k=2:iteration + + cell = glob.numberOfLayers(y,x,k); + while cell > 0 + + xco = [double(x-0.5)*glob.dx, double(x-0.5)*glob.dx, double(x+0.5)*glob.dx, double(x+0.5)*glob.dx]; + + allThick = sum(glob.thickness{y,x,k}(1:cell)); + oneThick = sum(glob.thickness{y,x,k}(1:cell-1)); + top = glob.strata(y,x,k-1)+allThick; + bottom = glob.strata(y,x,k-1)+oneThick; + fCode = glob.facies{y,x,k}(cell); % Note this is zero for no depositon, 9 for subaerial hiatus + + % Draw the in-situ production facies first + if fCode > 0 + zco = [bottom, top, top, bottom]; + faciesCol = [glob.faciesColours(fCode,1) glob.faciesColours(fCode,2) glob.faciesColours(fCode,3)]; + patch(xco, yco, zco, faciesCol,'EdgeColor','none'); + end + + cell = cell - 1; + end + end + end + + % Loop through iterations and draw timelines + for i=1:glob.timeLineCount + + k = glob.timeLineAge(i); + + if k <= iteration + for x = 1:long-1%glob.ySize-1 + % Draw a marker line across the top and down/up the side of a particular grid cell + xco = [x*glob.dx, (x+1)*glob.dx, (x+1)*glob.dx]; + zco = [glob.strata(y,x,k), glob.strata(y,x,k), glob.strata(y,x+1,k)]; + line(xco, [glob.ySize*glob.dx,glob.ySize*glob.dx,glob.ySize*glob.dx], zco, 'LineWidth',2, 'color', 'black'); + end + line([(x+1)*glob.dx,(x+2)*glob.dx], [double(glob.ySize)*glob.dx,double(glob.ySize)*glob.dx], [glob.strata(y,x+1,k), glob.strata(y,x+1,k)],'LineWidth',2, 'color', 'black'); + end + end +end \ No newline at end of file diff --git a/plotDipSections.m b/plotDipSections.m new file mode 100644 index 0000000..73c1216 --- /dev/null +++ b/plotDipSections.m @@ -0,0 +1,27 @@ +function plotDipSections(glob, stats, iteration) +% CHRONOSTRAT AND CROSS SECTION PLOT + + % ScreenSize is a four-element vector: [left, bottom, width, height]: + scrsz = get(0,'ScreenSize'); % vector + % position requires left bottom width height values. screensize vector is in this format 1=left 2=bottom 3=width 4=height + figure('Visible','on','Position',[20 scrsz(2)+20 scrsz(3)*0.9 scrsz(4)*0.6]); + + % plot includes cross section and chronostrat plotted on the same x axis + + xPos = round(glob.xSize / 2); % So line of section half-way across map grid + + waveProfilePlot = subplot('Position',[0.03, 0.90, 0.7, 0.09]); % position coords are x,y,width,height + plotWaveEnergyProfile(xPos, waveProfilePlot, glob); + + crossSectionPlot = subplot('Position',[0.03, 0.50, 0.7, 0.40]); + plotCrossSectionDip2D(xPos, crossSectionPlot, glob, iteration); + + chronoPlot = subplot('Position',[0.03, 0.075, 0.7, 0.4]); + plotChronostratSectionDip(xPos,chronoPlot, glob,iteration); + + SLPlot = subplot('Position',[0.75, 0.075, 0.1, 0.4]); + plotEustaticCurve(SLPlot,glob, iteration); + + PPlot = subplot('Position',[0.87, 0.075, 0.1, 0.4]); + plotFaciesVolumeCurves(PPlot, glob, stats, iteration); +end \ No newline at end of file diff --git a/plotEustaticCurve.m b/plotEustaticCurve.m new file mode 100644 index 0000000..3158398 --- /dev/null +++ b/plotEustaticCurve.m @@ -0,0 +1,23 @@ +function plotEustaticCurve(SLPlot, glob, iteration) + + subplot(SLPlot); + + % Sealevel curve + + % Force four ticks on the y axis and plot the y axis on the right hand side of the plot + set(SLPlot,'YAxisLocation','right'); + set(SLPlot,'YTick',[0 (glob.deltaT * iteration * 0.25) (glob.deltaT * iteration * 0.5) (glob.deltaT * iteration * 0.75) (glob.deltaT * iteration)]); + + for k=1:iteration-1 + % now plot the sea-level curve line for the same time interval + lineColor = [0.0 0.2 1.0]; + x = [double(glob.ySize)+glob.SL(k) double(glob.ySize)+glob.SL(k+1)]; + y = [double(k)*glob.deltaT double(k+1)*glob.deltaT]; + line(x,y, 'color', lineColor); + end + + xlabel('Eustatic Sealevel (m)','FontSize',11); + set(SLPlot,'YTick',[0 glob.deltaT * iteration]); + set(gca,'FontSize',11) + +end \ No newline at end of file diff --git a/plotFaciesVolumeCurves.m b/plotFaciesVolumeCurves.m new file mode 100644 index 0000000..6e3e9d9 --- /dev/null +++ b/plotFaciesVolumeCurves.m @@ -0,0 +1,16 @@ +function plotFaciesVolumeCurves(PPlot,glob, stats, iteration) +% Plot the production rate supply curve through time for each facies + + subplot(PPlot); + t = 1:iteration; + + for f = 1:glob.maxProdFacies * 2 % So plot the producing facies, and each transported facies derivative + xco = stats.totalFaciesVolume(t, f); + lineCol = [glob.faciesColours(f,1) glob.faciesColours(f,2) glob.faciesColours(f,3)]; + line(xco, t, 'color', lineCol, 'LineWidth', 2); + end + + xlabel('Facies volumes (m3)','FontSize',11); + set(gca,'FontSize',11); + grid on; +end \ No newline at end of file diff --git a/plotFinalGraphics.m b/plotFinalGraphics.m new file mode 100644 index 0000000..40aef1a --- /dev/null +++ b/plotFinalGraphics.m @@ -0,0 +1,13 @@ +function plotFinalGraphics(glob, stats, iteration) + + platMarginTraj = calculatePlatformMarginTrajectory(glob, iteration, round(glob.xSize / 2)); + + fprintf('Drawing the dip section and the chronostrat diagram...'); + plotDipSections(glob, stats, iteration); %cross and chrono + fprintf('Done\n'); + + fprintf('Drawing the 3D block diagram...'); + plot3DBlock(glob, iteration); %cross and chrono + fprintf('Done\n'); +end + diff --git a/plotInitialGraphics.m b/plotInitialGraphics.m new file mode 100644 index 0000000..fe8af69 --- /dev/null +++ b/plotInitialGraphics.m @@ -0,0 +1,135 @@ +function plotInitialGraphics(glob, graph, iteration) + + %cla % Clear all existing plots in the window + figure(graph.main); + + % Initial Subsidence map, position [left, bottom, width, height] all in range 0.0 to 1.0. + subsPlot = subplot('Position',[0.04 0.2 0.20 0.5]); + cla + reset(subsPlot); + p=surf(glob.subRateMap); + set(p,'LineStyle','none'); + xlabel('X Distance (grid points)'); + ylabel('Y Distance (grid points)'); + view([135, 45]); + colormap gray; + title('Subsidence map'); + + %% Initial bathymetry map + % position [left, bottom, width, height] all in range 0.0 to 1.0. + wdMapPlot = subplot('Position',[0.52 0.2 0.20 0.5]); + cla + reset(wdMapPlot); + p=surface(double(-glob.wd(:,:,1))); + set(p,'LineStyle','none'); + view([-85 25]); + grid on; + xlabel('X Distance (grid points)'); + ylabel('Y Distance (grid points)'); + view([135, 45]); + colormap gray; + title('Initial bathymetry'); + + %% Initial Facies map + % position [left, bottom, width, height] all in range 0.0 to 1.0. + faciesMapPlot = subplot('Position',[0.28 0.2 0.20 0.5]); + cla + reset(faciesMapPlot); + faciesM=zeros(glob.ySize,glob.xSize); + testX = 1:glob.xSize; + testY = 1:glob.ySize; + for y=1:glob.ySize + for x=1:glob.xSize + faciesM(y,x)=glob.facies{y,x,iteration}(1); + end + end + p=pcolor(testX, testY, double(faciesM)); + set(p,'LineStyle','none'); + view([0 90]); + xlabel('X Distance (grid points)'); + ylabel('Y Distance (grid points)'); + axis square + view([135, 45]); + colormap(glob.faciesColours); + title('Initial facies distribution'); + + %% Initial sea-level and carbonate production curves + wlPlot = subplot('Position',[0.04 0.8 0.6 0.15]); + cla + reset(wlPlot); + yyaxis left; + p = plot((1:glob.totalIterations) .* glob.deltaT, glob.SL(1:glob.totalIterations), 'linewidth', 3, 'color', [0, 0, 1]); + hold on + yyaxis right + p = plot((1:glob.totalIterations) .* glob.deltaT, glob.carbProdCurve(1:glob.totalIterations), 'linewidth', 3, 'color', [1, 0.1, 0]); + xlabel('Time (My)'); + ylabel('Carb prod modifier'); + yyaxis left; + ylabel('Eustatic Sea Level (m)'); + axis tight + + %% Production-depth curve plot + % initialise production depth curve plot + depthProdPlot = subplot('Position',[0.76 0.2 0.20 0.55]); + cla + reset(depthProdPlot); + + set(depthProdPlot, 'YDir', 'reverse'); + ylabel('Water depth (m)'); + xlabel('Production rates (m/My)'); + minMax = max(glob.prodRate/glob.deltaT); + axis([-(minMax*0.05) minMax*1.05 0 100]); + grid on; + depth=0:100; + + for j=1:glob.maxProdFacies + if glob.surfaceLight(j)>500 + %use this if you want to use Schlager's production profile + prodDepth = (glob.prodRate(j)/glob.deltaT) * tanh((glob.surfaceLight(j) * exp(-glob.extinctionCoeff(j) * depth))/ glob.saturatingLight(j)); + else + %use this if you want to use new production profile + prodDepth = (glob.prodRate(j)./glob.deltaT).* (1./ (1+((depth-glob.profCentre(j))./glob.profWidth(j)).^(2.*glob.profSlope(j))) ); + end + lineCol = [glob.faciesColours(j,1) glob.faciesColours(j,2) glob.faciesColours(j,3)]; + line(prodDepth, depth, 'color', lineCol, 'LineWidth',4); + end + + %% CA Rules + CARPlot=subplot('Position',[0.76 0.8 0.20 0.15]); + cla + reset(CARPlot); + radius=glob.CARules(1:glob.maxProdFacies,1); + bigRadius=max(radius(:)); + + % glob.faciesColours(1,2:4)=[27 170 185]./255; + % glob.faciesColours(2,2:4)=[197 202 32]./255; + % glob.faciesColours(3,2:4)=[140 84 161]./255; + + baseF=+0.5; + for f=1:glob.maxProdFacies + minS=glob.CARules(f,2); + maxS=glob.CARules(f,3); + minT=glob.CARules(f,4); + maxT=glob.CARules(f,5); + Col = [glob.faciesColours(f,1) glob.faciesColours(f,2) glob.faciesColours(f,3)]; + patch([minS maxS maxS minS],[baseF+0+0.1 baseF+0+0.1 baseF+0.5 baseF+0.5],Col,'linewidth',3); + patch([minT maxT maxT minT],[baseF+0.5 baseF+0.5 baseF+1-0.1 baseF+1-0.1],Col,'edgeColor','k','linewidth',3); + + posX=(minS+(maxS-minS)/2); + posY=baseF+.3; + text(posX,posY,'S'); + posX=(minT+(maxT-minT)/2); + posY=baseF+0.7; + text(posX,posY,'T'); + baseF=baseF+1; + end + axis([0 (((bigRadius*2+1))^2)-1 0.5 glob.maxProdFacies+.5]); + xTicks=[0:2:(((bigRadius*2+1))^2)-1]; + yTicks=[1:1:glob.maxProdFacies]; + + set(gca,'XTick',xTicks) + set(gca,'YTick',yTicks) + + xlabel('CA rules'); + ylabel('Facies'); +end diff --git a/plotWaveEnergyProfile.m b/plotWaveEnergyProfile.m new file mode 100644 index 0000000..9450e80 --- /dev/null +++ b/plotWaveEnergyProfile.m @@ -0,0 +1,13 @@ +function plotWaveEnergyProfile(xPos, waveProfilePlot, glob) + + subplot(waveProfilePlot); + + if strcmp(glob.waveRoutine,'on') == 1 + y = 1:glob.ySize; + waveEnergyProfile = glob.waveEnergy(y,xPos); + line(y, waveEnergyProfile, 'Color', [0.7,0.0,1.0], 'LineWidth',2, 'LineStyle','-.'); + end + + ylabel('Wave energy'); + grid on; +end \ No newline at end of file diff --git a/recordFaciesVolumes.m b/recordFaciesVolumes.m new file mode 100644 index 0000000..706ae5a --- /dev/null +++ b/recordFaciesVolumes.m @@ -0,0 +1,26 @@ +function stats = recordFaciesVolumes(glob, stats, iteration) + + faciesThicknesses = zeros(glob.ySize, glob.xSize, glob.maxFacies); + + % Loop across the model grid + for x = 1:glob.xSize + for y = 1:glob.ySize + for f = 1:glob.maxFacies % Loop through each facies + + % for the current iteration, each xy point for each facies f, sum the thicknesses of layers with facies code f + faciesThicknesses(y,x,f) = sum(glob.thickness{y,x,iteration}(glob.facies{y,x,iteration}==f)); + end + end + end + + % For each facies, find the sum thickness across the whole bvolume and + % multiply by grid cell size to get volume in m3 + gridCellArea = glob.dx * glob.dx; + f = 1:glob.maxFacies; + stats.totalFaciesVolume(iteration, f) = sum(sum(faciesThicknesses(:,:,f))) * gridCellArea; + + fprintf('Facies vols '); + for f=1:glob.maxProdFacies * 2; % all the produced facies and the transported facies + fprintf('%d:%4.3E ',f, stats.totalFaciesVolume(iteration, f)); + end +end \ No newline at end of file diff --git a/runCAModelGUI.m b/runCAModelGUI.m new file mode 100644 index 0000000..f0c326c --- /dev/null +++ b/runCAModelGUI.m @@ -0,0 +1,92 @@ +function [glob,stats,graph] = runCAModelGUI(glob, stats, graph) +% Run the model by executing all the functions specified in the processes input file + + order = 1; % rotates the order of facies neighbour checking in calculateFaciesCA to avoid bias for any one facies + iteration = 2; + + % initialize movie + %[graph,theMovie]=initializeMovie(graph); + + while iteration <= glob.totalIterations + + fprintf('It:%d EMT %4.3f ', iteration, glob.deltaT * iteration); + + glob.strata(:,:,iteration) = glob.strata(:,:,iteration-1); + + glob = calculateSubsidence(glob, iteration); + + glob = calculateWaterDepth(glob, iteration); % NB water depth required for facies distrib so needs to be calculated first + + if strcmp(glob.waveRoutine,'on') == 1 + glob = calculateWaveEnergy(glob,iteration); + end + + % Choose if the carbonate concentration routine should be included------------------------------------- + if strcmp(glob.concentrationRoutine,'on') == 1 + glob = diffuseConcentration(glob,iteration); + end + + % Calculate facies distribtion whichever CA routine has been selected ---------------------------------------------------------------------- + switch glob.CARoutine + case 'simpleWave' + glob = calculateFaciesWaveDependentSimple(glob, iteration); + case 'orderedCA' + glob = calculateFaciesCARotationOrder(glob, iteration, order); + case 'productionCA' + glob = calculateFaciesCAProdDependent(glob, iteration); + case 'waveCA' + glob = calculateFaciesCAWaveDependent(glob, iteration); + case 'waveHybridCA' + glob = calculateFaciesCAWaveDependentHybrid(glob, iteration); + otherwise + glob = calculateFaciesCARotationOrder(glob, iteration, order); % Simplest CA option is the default + end + + % input and diffusionally transport siliciclastics--------------------------------------------------------------------- + if strcmp(glob.siliciclasticsRoutine,'on') == 1 + glob = diffuseSiliciclastic (glob,iteration); + end + + % Carbonate production----------------------------------------------------------------------------------- + glob = calculateProduction(glob, iteration); + + % Choose carbonate transportation-------------------------------------------------------------------- + if strcmp(glob.transportationRoutine,'steepestSlope') == 1 + glob = calculateTransport(glob, iteration); + else + glob = calculateTransportCrossPlatform (glob, iteration); + end + + % Choose if the soil deposition routine should be included--------------------------------------------------------------- + checkProcess=strcmp(glob.soilRoutine,'on'); + if checkProcess==1 + glob = calculateSoilDeposition (glob,iteration); + end + + stats = recordFaciesVolumes(glob, stats, iteration); + + %CarboCAT the movie------------------------------------------------------------------------------------- + %[graph,theMovie]=recordMovie(glob,iteration,graph,theMovie,subs); + %[graph,theMovie]=recordMovieSimple(glob,iteration,graph,theMovie,subs); + + % Control cycle through facies for neighbour checking + order = order + 1; + if order > glob.maxProdFacies; order = 1; end + iteration = iteration + 1; + fprintf('\n'); + end + + %close (theMovie); + + % Reverse the effect of the final increment to avoid problems with array overflows in graphics etc + if iteration > glob.totalIterations + iteration = iteration - 1; + end + + fprintf('Model complete after %d iterations and ready to plot\n',iteration); + + % movie(mapMovie,1,1); + % movie2avi + %graph = finalGraphics(glob, stats, graph, iteration); + % graph = movie(glob, stats, graph, iteration); +end