forked from DomHu/iTURBO2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
iturbo2script_multipleSims_dissolution.m
201 lines (174 loc) · 8.07 KB
/
iturbo2script_multipleSims_dissolution.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
function iturbo2script_multipleSims_dissolution(data, data50, data90, data100, nocarriers, noexps, expname, TM, plot_min)
%% MATLAB script to run TURBO2 experiment with input data
%% RUN three different dissolution experiments AND PLOT results in same Fig.
% data = matrix of required data (age, mxl, abu, iso) without dissoution (abundance change)
% data50 = matrix of required data with 50% dissoution (50% abundance change)
% data90 = matrix of required data with 90% dissoution (90% abundance change)
% data100 = matrix of required data with 100% dissoution (100% abundance change)
% nocarriers = number of carriers to be measured
% noexps = number of experiments for SA
% expname = experiment name
% TM = true or false for using trnasition matrix or homogeneous mixing
% plot_min = true/false; just for plotting purposes -- to make sure the min point is plotted
% age = age of sediment layer down core
% mxl = series of mixed layer thicknesses (zbio) down core
% abu = series of abundances of carrier type 1 down core
% iso = original isotope signature of both carrier types 1 and 2
c = @cmu.colors; % shortcut function handle
age = data(:,1); % here asme for all data-files
mxl1 = data(:,2); % here asme for all data-files
abu = data(:,3);
abu50 = data50(:,3);
abu90 = data90(:,3);
abu100 = data100(:,3);
iso = data(:,4); % here asme for all data-files
lngth = length(data(:,1));
numb = nocarriers; % number of carriers to be measured
exps = noexps; % number of different experiments
%%
for i = 1:exps
[oriabu(i,:,:),bioabu(i,:,:),oriiso(i,:,:),bioiso(i,:,:)] = iturbo2_plus_TM(abu,iso,mxl1,numb, TM);
[oriabu2(i,:,:),bioabu2(i,:,:),oriiso2(i,:,:),bioiso2(i,:,:)] = iturbo2_plus_TM(abu50,iso,mxl1,numb, TM);
[oriabu3(i,:,:),bioabu3(i,:,:),oriiso3(i,:,:),bioiso3(i,:,:)] = iturbo2_plus_TM(abu90,iso,mxl1,numb, TM);
[oriabu4(i,:,:),bioabu4(i,:,:),oriiso4(i,:,:),bioiso4(i,:,:)] = iturbo2_plus_TM(abu100,iso,mxl1,numb, TM);
end
%%
% variable for mean results of no dissolution
mean_bioabu1_diss0 = zeros(1,lngth);
mean_bioabu2_diss0 = zeros(1,lngth);
mean_bioiso1_diss0 = zeros(1,lngth);
mean_bioiso2_diss0 = zeros(1,lngth);
% variable for mean results of 50% dissolution
mean_bioabu1_mxl2 = zeros(1,lngth);
mean_bioabu2_mxl2 = zeros(1,lngth);
mean_bioiso1_mxl2 = zeros(1,lngth);
mean_bioiso2_mxl2 = zeros(1,lngth);
% variable for mean results of 90% dissolution
mean_bioabu1_mxl3 = zeros(1,lngth);
mean_bioabu2_mxl3 = zeros(1,lngth);
mean_bioiso1_mxl3 = zeros(1,lngth);
mean_bioiso2_mxl3 = zeros(1,lngth);
% variable for mean results of 100% dissolution
mean_bioabu1_diss100 = zeros(1,lngth);
mean_bioabu2_diss100 = zeros(1,lngth);
mean_bioiso1_diss100 = zeros(1,lngth);
mean_bioiso2_diss100 = zeros(1,lngth);
%%
mxltext = num2str(mean(mxl1));
numbtxt = num2str(numb,2);
expstxt = num2str(exps,2);
abutxt = num2str(max(abu),2);
set(0,'DefaultAxesFontSize',16)
%% Plot isotope plots only and just for species 1 (no abundance change here)
fig1 = figure;
hold on
for i = 1:exps
% dissolution 90%
plot(1:lngth,bioiso3(i,:,1), 'Color', [0.8 0.8 1.0],'Linewidth',1.5)
mean_bioiso1_mxl3 = mean_bioiso1_mxl3+bioiso3(i,:,1);
% dissolution 50%
plot(1:lngth,bioiso2(i,:,1), 'Color', [0.8 1.0 0.8],'Linewidth',1.5)
mean_bioiso1_mxl2 = mean_bioiso1_mxl2+bioiso2(i,:,1);
% dissolution 0%
plot(1:lngth,bioiso(i,:,1), 'Color', [1.0 0.8 0.8],'Linewidth',1.5)
mean_bioiso1_diss0 = mean_bioiso1_diss0+bioiso(i,:,1);
% dissolution 100%
plot(1:lngth,bioiso4(i,:,1), 'Color',c('carrot orange'),'Linewidth',1.5)
mean_bioiso1_diss100 = mean_bioiso1_diss100+bioiso4(i,:,1);
end
plot(1:lngth,oriiso(1,:,1),'k','Linewidth',2.0) % plot original iso for species 1
% mxl3
mean_bioiso1_mxl3 = mean_bioiso1_mxl3/exps;
plot(1:lngth,mean_bioiso1_mxl3, '-b','Linewidth',2.0)
% mxl2
mean_bioiso1_mxl2 = mean_bioiso1_mxl2/exps;
plot(1:lngth,mean_bioiso1_mxl2, '-g','Linewidth',2.0)
% mxl1
mean_bioiso1_diss0 = mean_bioiso1_diss0/exps;
plot(1:lngth,mean_bioiso1_diss0, '-r','Linewidth',2.0)
% dissolution 100%
mean_bioiso1_diss100 = mean_bioiso1_diss100/exps;
plot(1:lngth,mean_bioiso1_diss100, '-','Color',c('deep carrot orange'),'Linewidth',2.0)
set(gca,'YDir','Reverse','XGrid','On','YGrid','On','Box','On', 'XLim',[0,200], 'YLim',[1.0,3.0],'YTick',[1.0 1.5 2.0 2.5 3.0])
xlabel('Core depth (cm) ');
ylabel('\delta^{18}O');
printfilename = [expname,'_',abutxt,'abu_',numbtxt,'carriers_',expstxt,'Exps'];
% check if output directory exists -- if not create it:
if ~(exist('output/mat','dir') == 7), mkdir('output/mat'); end
save(['output/mat/',printfilename,'.mat'],'printfilename', 'lngth','bioiso','bioiso2','bioiso3','bioiso4', 'oriiso', 'mean_bioiso1_diss0', 'mean_bioiso1_mxl2', 'mean_bioiso1_mxl3', 'mean_bioiso1_diss100','expname', 'exps', 'bioabu','bioabu2','bioabu3','bioabu4', 'oriabu', 'oriabu2', 'oriabu3', 'oriabu4', 'mean_bioabu1_diss0', 'mean_bioabu1_mxl2', 'mean_bioabu1_mxl3', 'mean_bioabu1_diss100')
print(fig1, '-depsc', ['output/',printfilename]); % save figure in extra output folder
%% Plot abundance of species 1
fig2 = figure;
hold on
for i = 1:exps
% 0% dissolution
plot(1:lngth,bioabu(i,:,1), 'Color', [1.0 0.8 0.8],'Linewidth',1.5)
mean_bioabu1_diss0 = mean_bioabu1_diss0+bioabu(i,:,1);
% 50% dissolution
plot(1:lngth,bioabu2(i,:,1), 'Color', [0.8 1.0 0.8],'Linewidth',1.5)
mean_bioabu1_mxl2 = mean_bioabu1_mxl2+bioabu2(i,:,1);
% 90% dissolution
plot(1:lngth,bioabu3(i,:,1), 'Color', [0.8 0.8 1.0],'Linewidth',1.5)
mean_bioabu1_mxl3 = mean_bioabu1_mxl3+bioabu3(i,:,1);
% 100% dissolution
plot(1:lngth,bioabu4(i,:,1), 'Color',c('carrot orange'),'Linewidth',1.5)
mean_bioabu1_diss100 = mean_bioabu1_diss100+bioabu4(i,:,1);
end
% 100% dissolution
mean_bioabu1_diss100 = mean_bioabu1_diss100/exps;
plot(1:lngth,mean_bioabu1_diss100, 'Color',c('deep carrot orange'),'Linewidth',2.0)
% 90% dissolution
mean_bioabu1_mxl3 = mean_bioabu1_mxl3/exps;
plot(1:lngth,mean_bioabu1_mxl3, '-b','Linewidth',2.0)
% 50% dissolution
mean_bioabu1_mxl2 = mean_bioabu1_mxl2/exps;
plot(1:lngth,mean_bioabu1_mxl2, '-g','Linewidth',2.0)
% 0% dissolution
mean_bioabu1_diss0 = mean_bioabu1_diss0/exps;
plot(1:lngth,mean_bioabu1_diss0, '-r','Linewidth',2.0)
% plot one of the original abu 0% dissolution
line_fewer_markers(1:lngth,oriabu(1,:,1), 20,'--ko',...
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','r',...
'MarkerSize',6);
% plot one of the original abu 50% dissolution
line_fewer_markers(1:lngth,oriabu2(1,:,1),20,'--ko',...
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',6);
% plot one of the original abu 90% dissolution
line_fewer_markers(1:lngth,oriabu3(1,:,1), 20,'--ko',...
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','b',...
'MarkerSize',6);
% plot one of the original abu 100% dissolution
line_fewer_markers(1:lngth,oriabu4(1,:,1), 20,'--ko',...
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor',c('deep carrot orange'),...
'MarkerSize',6);
% plot minimum number of particles
if(plot_min)
% find min
idmin50perc = find(oriabu2(1,:,1) == min(oriabu2(1,:,1)));
idmin90perc = find(oriabu3(1,:,1) == min(oriabu3(1,:,1)));
idmin100perc = find(oriabu4(1,:,1) == min(oriabu4(1,:,1)));
plot(idmin50perc, oriabu2(1,idmin50perc,1),'o-', 'LineWidth',2, 'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',6)
plot(idmin90perc, oriabu3(1,idmin90perc,1),'o-', 'LineWidth',2, 'MarkerEdgeColor','k',...
'MarkerFaceColor','b',...
'MarkerSize',6)
plot(idmin100perc, oriabu4(1,idmin100perc,1),'o-', 'LineWidth',2, 'MarkerEdgeColor','k',...
'MarkerFaceColor',c('deep carrot orange'),...
'MarkerSize',6)
end
set(gca,'XGrid','On','YGrid','On', 'XLim',[0,200], 'YLim',[0,1000], 'YTick',[0 200 400 600 800 1000])
xlabel('Core depth (cm) ');
ylabel('Number of Particles');
title('Abundance of Species 1')
printfilename = [expname,'_',abutxt,'abu_',numbtxt,'carriers_',expstxt,'Exps_ABU'];
print(fig2, '-depsc', ['output/',printfilename]); % save figure in extra output folder