% This script does analyses of GISP2 bedrock core Be-10 and Al-26 data.
% Written by Greg Balco
% Berkeley Geochronology Center
% May-August 2016
%
% Note that this requires the MATLAB optimization toolbox. 
%
% It won't run in Octave due to extensive plotting calls. I am using MATLAB
% R2015b. I imagine it will run in some earlier versions.  
% 
% I suggest running this section-by-section, because only one figure is 
% allowed to be alive at a time. Thus, if you run the whole thing you will
% not see the results. 
%
% Also note that most of the sections are
% wrapped with an if 1...end block so they can be turned on and off. 
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License, version 2,
% as published by the Free Software Foundation (www.fsf.org).

% Clean up

clear all; close all;  

%% Enter data

% A total of 8 analyses were made, which overlap in depth somewhat.
% Also, the quartz yield was variable from the different core segments that
% were mixed together. 
%
% Here organize all these data in a sensible fashion. 

% These data updated May 2016. 

% Be-10 analysis on 2/3
s{1}.name = 'GISP23';
s{1}.nuclide = 10; % 10 = Be-10; 26 = Al-26
s{1}.td = [0.22 0.65]; % Top depths of core segments (m)
s{1}.bd = [0.315 0.8]; % Bottom depths of core segments (m)
s{1}.mq = [3.6127 9.0427]; % g of quartz from each core segment
s{1}.N = 2.480e4; % Measured nuclide concentration
s{1}.dN = 7.400e2; % Uncertainty therein

% Be-10 analysis on 4
s{2}.name = 'GISP4';
s{2}.nuclide = 10;
s{2}.td = 0.8;
s{2}.bd = 0.994;
s{2}.mq = 17.8793;
s{2}.N = 1.797e4; 
s{2}.dN = 7.20e2;

% Be-10 analysis on 5/6
s{3}.name = 'GISP56';
s{3}.nuclide = 10;
s{3}.td = [1.02 1.10];
s{3}.bd = [1.10 1.135];
s{3}.mq = [6.4583 2.8823];
s{3}.N = 1.460e4; 
s{3}.dN = 1.00e3;

% Be-10 analysis on 7
s{4}.name = 'GISP7';
s{4}.nuclide = 10;
s{4}.td = 1.14;
s{4}.bd = 1.285;
s{4}.mq = 8.8088;
s{4}.N = 1.180e4; 
s{4}.dN = 8.00e2;

% Be-10 analysis on 8/9
s{5}.name = 'GISP89';
s{5}.nuclide = 10;
s{5}.td = [1.285 1.42];
s{5}.bd = [1.40 1.55];
s{5}.mq = [5.0505 3.8182];
s{5}.N = 9.800e3; 
s{5}.dN = 4.90e2;

% Al-26 analysis on 2/3/4
s{6}.name = 'GISP234';
s{6}.nuclide = 26;
s{6}.td = [0.22 0.65 0.8];
s{6}.bd = [0.315 0.8 0.994];
s{6}.mq = [3.6127 9.0437 17.8793];
s{6}.N = 8.80e4; 
s{6}.dN = 5.20e3;

% Al-26 analysis on 5/6/7
s{7}.name = 'GISP567';
s{7}.nuclide = 26;
s{7}.td = [1.02 1.10 1.14];
s{7}.bd = [1.10 1.135 1.285];
s{7}.mq = [6.4583 2.8823 8.8088];
s{7}.N = 5.49e4; 
s{7}.dN = 5.38e3;

% Al-26 analysis on 8/9 -- this is preliminary measurement whose Al yield
% could not be verified; we don't use this in the model fitting, but you
% can uncomment it to plot it on the figure if you want. 
% s{8}.name = 'GISP89';
% s{8}.nuclide = 26;
% s{8}.td = [1.285 1.42];
% s{8}.bd = [1.40 1.55];
% s{8}.mq = [5.0505 3.8182];
% s{8}.N = 3.04e4; 
% s{8}.dN = 4.74e3;


%% Plot data

% This generates a plot showing all the measurements that were made. 

ff = plotGISPdata(s);

% Trim extra stuff for appearance purposes 
aa = get(gcf,'children');
delete(aa(1));
set(aa(2),'pos',[0.1300    0.1100    0.7750    0.8150]);
grid on;


%% Define constants and other inputs necessary for forward model calculations. 

% Define constants. 

consts.R2610q = 6.75; % Spallogenic 26/10 ratio in quartz
consts.P10q_St = 4.1; % This is NENA (Balco et al., 2009) after accounting
% for reduced muon production.
% Note that the uncertainty here is much less important than the 
% uncertainty in elevation and air pressure during past exposure periods 
% when the ice sheet was absent. 
consts.delP10q_St = 0.19; % Corresponding uncert. Not particularly relevant. 

% Muon interaction cross-sections

% Be-10, from John Stone - BCO core fit
consts.mc10q.Natoms = 2.006e22;
consts.mc10q.k_neg = 0.704 * 0.1828 * 0.0011;
consts.mc10q.delk_neg = consts.mc10q.k_neg.*0.1; % Estimated
consts.mc10q.sigma190 = 0.81e-30.*(190^0.75); % Stone reported sigma0
consts.mc10q.delsigma190 = consts.mc10q.sigma190.*0.1; % estimated

% Al-26 from John Stone - BCO core fit
consts.mc26q.Natoms = 1.003e22;
consts.mc26q.k_neg = 0.296 * 0.6559 * 0.0084;
consts.mc26q.delk_neg = consts.mc26q.k_neg.*0.1; % Estimated
consts.mc26q.sigma190 = 13.6e-30.*(190^0.75); % Stone reported sigma0
consts.mc26q.delsigma190 = consts.mc26q.sigma190.*0.1; % estimated

% Also define decay constants

consts.l10 = -log(0.5)./1.387e6; % Chmeleff/Korschinek value
dldt = -log(0.5).*(1.387e6^-2);
consts.dell10 = sqrt((dldt.*0.012e6)^2); 
consts.l26 = 9.83e-7; % Compatible with Nishiizumi (2004)
consts.dell26 = 2.5e-8;

% Convert linear depth in core to mass depth. 

rho = 2.7;
for a = 1:length(s);
    s{a}.tdg = s{a}.td.*100.*rho;
    s{a}.bdg = s{a}.bd.*100.*rho;
end;

% Latitude and elevation of core site
site_lat =  72.5796;
site_long = -38.4592;
site_elv = 0; % Assume zero for simplicity. Speculative. 
 site_elv = 200; % Uncomment to try a different elevation. 
% Atmospheric pressure at site
site_p = NCEPatm_2(site_lat,site_long,site_elv); % site air pressure
% Of course, that is even more speculative because the atmospheric 
% circulation would be totally different if the ice sheet wasn't there. 

% Build and load muon profile
% build_muon_profile.m builds a production rate profile defined on a grid
% for efficient integration later.
m = build_muon_profile(site_p,consts,1);

% Define production rate info
SFsp = stone2000(site_lat,site_p,1); % scaling factor

% Build a data structure with production rate information
p.P10sp = consts.P10q_St .* SFsp; % Be-10 spallation production rate at surface 
p.P26sp = p.P10sp.*consts.R2610q; % Al-26 spallation production rate at surface

% Attenuation
p.Lsp = 140; % g/cm2, appropriate for high latitude. 

% Set optimization rules... 
opts1 = optimset('fminsearch');
opts1 = optimset(opts1,'TolFun',1e-3,'TolX',1e-3);

opts2 = optimset('fmincon');
opts2 = optimset(opts2,'TolFun',1e-3,'TolX',1e-3,'display','iter');


%% Try to fit production profile without regard to Al/Be production ratio, 
% assuming surface is zero depth. This demonstrates that the apparent
% attenuation length of the measured concentrations is much longer than
% that for spallogenic production near the surface. 

if 1;

% Define L 
% Allow N10,0 and N26,0 to float

% This generates the first part of Extended Data Figure 1. 

delete(ff);
ff = plotGISPdata(s);
% Trim extra stuff
delete(ff(3));
set(ff(2),'pos',[0.1300    0.1100    0.7750    0.8150]);
grid on;

mask = [1:7]; 

[optX,fmin] = fminsearch(@(X) fitSpMu(X,s,p,m,mask,ff),[20 13],opts1)

disp(['Production profile fixed to zero depth. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-2))]);

% Clean up figure for publication
%set(findobj('parent',gca,'color','g'),'color','k');
%set(gca,'fontsize',12);
%title([]);

end;

%% Try to fit with spallation and muons fixed to production profile, but floating depth
% Al-26/Be-10 is also allowed to float without limits here. 
% This demonstrates that you can fit the profile with a single period of 
% exposure without any erosion, if there was additional mass covering the
% bedrock surface. 

if 1;

% This generates Figure S3.
close all;
ff = plotGISPdata(s);
% Trim extra stuff
delete(ff(3));
set(ff(2),'pos',[0.1300    0.1100    0.7750    0.8150]);
grid on;

mask = [1:7]; 

[optX,fmin] = fminsearch(@(X) fitSpMuZ(X,s,p,m,mask,ff),[150 80 350],opts1)

disp(['Production profile with floating depth. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-3))]);

% Clean up figure for publication
%set(findobj('parent',gca,'color','g'),'color','k');
%set(gca,'fontsize',12);
%title([]);

end;

%% Do a Monte Carlo simulation for formal uncertainty in cover depth 
% estimate. 

if 0;

close all;
ff = plotGISPdata(s);
numits = 100;
mask = 1:7;

for a = 1:numits;
    s1 = s;
    for b = 1:length(s);
        % Generate random realization of nuclide concentrations
        s1{b}.N = s{b}.N + s{b}.dN.*randn(1);
    end;
    [optX,fmin] = fminsearch(@(X) fitSpMuZ(X,s1,p,m,mask,ff),[150 80 350],opts1);
    thisz(a) = optX(3);
    thisfmin(a) = fmin;
    disp([int2str(a) ' -- ' sprintf('%0.1f',optX(3))]);
end;

disp(['Best-fit cover thickness is ' sprintf('%0.1f',mean(thisz)) ' +/- ' sprintf('%0.1f',std(thisz))]);


end;


%% Two-stage exposure history

% Now try to fit data with two-stage exposure history at floating depth.
% Three parameters: Exposure time, burial time, depth below actual surface

% This generates the upper half of Extended Data Figure 2. 

if 1;

close all;
ff = plotGISPdata(s);
    
mask = [1:7];

[optX,fmin] = fminsearch(@(X) fit_tE_tB_zz(X,s,p,m,mask,ff),[300 1 300],opts1)

disp(['Two-stage exposure-burial, no erosion, cover. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-3)) ' for ' int2str(length(mask)-3) ' DOF']);


% Also plot on splot

pred = fit_tE_tB_zz(optX,s,p,m,mask,[],1);
GISP2_splot(s,optX(3),m,p,consts,pred);


% Clean up figure for publication
set(findobj('color','g'),'color','k');
set(findobj('markerfacecolor','g'),'markerfacecolor','k');
set(get(gcf,'children'),'fontsize',12);

end;

%% Do a MCS to determine the uncertainty in the above. 

if 0; 

mask = [1:7];
numits = 100;
MCtexp = zeros(1,numits); MCtbur = MCtexp; MCcovr = MCtexp; fmins = MCtexp;
for a = 1:length(s);
    MCN(a,:) = randn(1,numits).*s{a}.dN + s{a}.N;
end;
for a = 1:numits;
    s1 = s;
    for b = 1:length(s1);
        s1{b}.N = MCN(b,a);
    end;
        x0 = rand(1,3).*[200 1 200] + [50 0.5 150];
        [optX,fmin] = fminsearch(@(X) fit_tE_tB_zz(X,s1,p,m,mask,[]),x0,opts1);
        MCtexp(a) = optX(1);
        MCtbur(a) = optX(2);
        MCcovr(a) = optX(3);
        fmins(a) = fmin;
        disp(['Completed ' int2str(a) ': tE = ' sprintf('%0.0f',optX(1)) '; tB = ' sprintf('%0.2f',optX(2))]);
end;

 
% linearly approximate external uncertainty 
ll = consts.l26 - consts.l10; dll = sqrt(consts.dell26.^2 + consts.dell10.^2);
RR0 = exp(-ll.*1e6.*mean(MCtbur)); 
dtdR = -1./(ll.*RR0);
dRR0 = 1e6.*std(MCtbur) ./ dtdR;
dtdll = log(RR0).*(ll.^-2);
deltext = sqrt( (dtdR.*dRR0).^2 + (dtdll.*dll).^2)
disp(['Burial age: ' sprintf('%0.3f',mean(MCtbur)) ' +/- ' sprintf('%0.3f',std(MCtbur)) ' (' sprintf('%0.3f',deltext./1e6) ') Ma']);


end;

%% Steady erosion

% Try to fit with steady erosion followed by burial. No sediment cover. 
% Two parameters: steady state erosion rate and burial time

% This generates Figures S7 and S8. 

if 1;
    
close all;
ff = plotGISPdata(s);

mask = [1:7];

[optX,fmin] = fminsearch(@(X) fit_steadyE_tB_no_z(X,s,p,m,mask,ff),[30 1],opts1)

disp(['Steady erosion-burial, no cover. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-2)) ' for ' int2str(length(mask)-2) ' DOF']);

% Also plot on splot

pred = fit_steadyE_tB_no_z(optX,s,p,m,mask,[],1);
GISP2_splot(s,0,m,p,consts,pred);

% Clean up figure for publication
set(findobj('color','g'),'color','k');
set(findobj('markerfacecolor','g'),'markerfacecolor','k');
set(get(gcf,'children'),'fontsize',12);

end;

%% Two-stage mixed exposure-erosion scenarios; parameter space random search
% This just identifies scenarios that are OK to show that there are lots of
% acceptable combinations. It takes a long time to run. 

if 0;
    
    % This generates Figure 9. It takes a while to run. 
    close all;
    ff = plotGISPdata(s);

    mask = [1:7];

    numits = 100;
    start_te = rand(1,numits).* 1000; % Ka
    start_tb = rand(1,numits) + 0.5; % Ma
    start_zz = rand(1,numits).* 400; % g/cm2
    start_ee = rand(1,numits).*34; % m/Myr

    optxs = zeros(numits,4);
    fmins = zeros(1,numits);

    figure(2); 
    subplot(2,2,1); xlabel('texp(ka)'); ylabel('Erosion rate (m/Myr)');
    subplot(2,2,2); xlabel('texp(ka)'); ylabel('Cover (g/cm2)');
    subplot(2,2,3); xlabel('Erosion rate (m/Myr)'); ylabel('Cover (g/cm2)');
    subplot(2,2,4); xlabel('Erosion rate (m/Myr)'); ylabel('Burial time (Ma)');

    opts3 = optimset(opts2,'display','none');

    for a = 1:numits;
        X0 = [start_te(a) start_tb(a) start_ee(a) start_zz(a)];
        [optxs(a,:), fmins(a)] = fmincon(@(X) fit_tE_tB_e_zz(X,s,p,m,mask,[]),X0,[],[],[],[],[0 0.5 0 0],[Inf 1.5 35 1000],[],opts3);
        figure(2); subplot(2,2,1); plot(optxs(a,1),optxs(a,3),'ro'); hold on;
        subplot(2,2,2); plot(optxs(a,1),optxs(a,4),'ro'); hold on;
        subplot(2,2,3); plot(optxs(a,3),optxs(a,4),'ro');hold on;
        subplot(2,2,4); plot(optxs(a,3),optxs(a,2),'ro'); hold on; drawnow;
    end;

    %% Nicer plotting code for this -- this generates Figure S9

    ok = find(fmins < 2.6*3); % Exclude poor fits
    figure;
    subplot(2,2,1); 
    plot(optxs(ok,1),optxs(ok,3),'bo','markerfacecolor',[0.5 0.5 1])
    xlabel('Exposure time(ka)'); ylabel('Erosion rate (m/Myr)');
    axis([200 800 0 35]); grid on;

    subplot(2,2,2);
    plot(optxs(ok,1),optxs(ok,4),'bo','markerfacecolor',[0.5 0.5 1]);
    xlabel('Exposure time(ka)'); ylabel('Cover (g/cm2)');
    axis([200 800 0 400]); grid on;

    subplot(2,2,3);
    plot(optxs(ok,3),optxs(ok,4),'bo','markerfacecolor',[0.5 0.5 1]);
    xlabel('Erosion rate (m/Myr)'); ylabel('Cover (g/cm2)');
    axis([0 35 0 400]); grid on;

    subplot(2,2,4);
    plot(optxs(ok,1),optxs(ok,2),'bo','markerfacecolor',[0.5 0.5 1]);
    xlabel('texp(ka)');ylabel('Burial time (Ma)');
    axis([200 800 0.5 1.5]); grid on;

end;

%% Now consider no-erosion, multiple exposure-burial scenarios. 

%% Dynamic glacial-interglacial steady state scenario. 
% Only free parameter here is length of interglacials; cycle assumed 100 ka
% This generates the first part of Extended Data Figure 4

if 1;
    
close all;
ff = plotGISPdata(s);

mask = [1:7];

[optX,fmin] = fminsearch(@(X) fit_dynamic_ss(X,s,p,m,mask,350,1e5,ff,0),10,opts1)

disp(['Dynamic glacial-interglacial steady state, 100 ka. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-1)) ' for ' int2str(length(mask)-1) ' DOF']);

pred = fit_dynamic_ss(optX,s,p,m,mask,350,1e5,ff,1);
GISP2_splot(s,350,m,p,consts,pred);

% Clean up figure for publication
set(findobj('color','g'),'color','k');
set(findobj('markerfacecolor','g'),'markerfacecolor','k');
figure(1);
set(get(gcf,'children'),'fontsize',12);
figure(2);
set(get(gcf,'children'),'fontsize',12);

end;

%% show 1-param fit graphically

if 0;
    
    close all;
    ff = plotGISPdata(s);
    mask = [1:7];
    xs = 1:30;
    for a = 1:length(xs);
        optX(a) = fit_dynamic_ss(xs(a),s,p,m,mask,350,1e5,ff,0);
    end;
    
    figure;
    semilogy(xs,optX);
    clear optX;
end;


%% Various glacial-interglacial scenarios. 

if 0;

mask = [1:7];

close all;
ff = plotGISPdata(s);

% Scenario with all late Pleistocene interglacials, length adjustable, and
% a period of exposure at 1.1 Ma. 

mids = fliplr([125 235 335 415 525 615 705 785 865 935 965 1005 1065]) .* 1000;

%[optX,fmin] = fminsearch(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[150 10],opts1)

opts2 = optimset('fmincon');
opts2 = optimset(opts2,'TolFun',1e-4,'TolX',1e-4,'display','iter');

[optX,fmin] = fmincon(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[200 5],[],[],[],[],[0 0],[Inf 80],[],opts2)

disp(['Mid-late Pleistocene scenario 0. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-2)) ' for ' int2str(length(mask)-2) ' DOF']);

% Of course, if you think about it for three seconds, you will realize that
% that scenario just goes straight to zero interglacial length and 
% converges on Model 1. So that exercise was not productive. Ignore. 

end;
    
%% Scenario with all late Pleistocene interglacials and some 41 kyr cycles 
% tacked on at the beginning.

% This generates the second part of Extended Data Figure 4

if 1;
    
close all;
ff = plotGISPdata(s);


mask = [1:7];

mids = fliplr([125 235 335 415 525 615 705 785 865 935 965 1005 1065]) .* 1000;

adds = mids(1) + fliplr(cumsum(zeros(1,10)+41000));
mids = [adds mids];

opts2 = optimset('fmincon');
opts2 = optimset(opts2,'TolFun',1e-4,'TolX',1e-4,'display','iter');

[optX,fmin] = fmincon(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[200 5],[],[],[],[],[0 0],[Inf 80],[],opts2)

disp(['Mid-late Pleistocene scenario 1. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-2)) ' for ' int2str(length(mask)-2) ' DOF']);


% Clean up figure for publication
set(findobj('color','g'),'color','k');
set(findobj('markerfacecolor','g'),'markerfacecolor','k');
figure(1);
set(get(gcf,'children'),'fontsize',12);

end;

%% Scenario with largest Pleistocene interglacials and some 41 kyr cycles 
% tacked on at the beginning. 

if 1;
close all;
ff = plotGISPdata(s);

% This generates the third part of Extended Data Figure 4

mask = [1:7];

mids = fliplr([125 335 415 1065]) .* 1000;% only 9/11/31

adds = mids(1) + fliplr(cumsum(zeros(1,10)+41000));
mids = [adds mids];

%[optX,fmin] = fminsearch(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[150 10],opts1)

opts2 = optimset('fmincon');
opts2 = optimset(opts2,'TolFun',1e-4,'TolX',1e-4,'display','iter');

[optX,fmin] = fmincon(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[200 5],[],[],[],[],[0 0],[Inf 80],[],opts2)

disp(['Mid-late Pleistocene scenario 2. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-2)) ' for ' int2str(length(mask)-2) ' DOF']);

% Clean up figure for publication
set(findobj('color','g'),'color','k');
set(findobj('markerfacecolor','g'),'markerfacecolor','k');
figure(1);
set(get(gcf,'children'),'fontsize',12);

end;

%% Scenario with largest Pleistocene interglacials and a bunch of 164 kyr cycles 
% tacked on at the beginning, i.e. ice sheet disappears in 1 of 4 41 ka cycles. 

if 1;
close all;
ff = plotGISPdata(s);

% This generates the last part of Extended Data Figure 4.

mask = [1:7];

mids = fliplr([125 335 415 1065]) .* 1000;% only 9/11/31

adds = mids(1) + fliplr(cumsum(zeros(1,10)+164000));
mids = [adds mids];

%[optX,fmin] = fminsearch(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[150 10],opts1)

opts2 = optimset('fmincon');
opts2 = optimset(opts2,'TolFun',1e-4,'TolX',1e-4,'display','iter');

[optX,fmin] = fmincon(@(X) fit_midp1(X,s,p,m,mask,mids,ff),[200 5],[],[],[],[],[0 0],[Inf 80],[],opts2)

disp(['Mid-late Pleistocene scenario 3. Reduced chi-squared is ' sprintf('%0.2f',fmin./(length(mask)-2)) ' for ' int2str(length(mask)-2) ' DOF']);

% Clean up figure for publication
set(findobj('color','g'),'color','k');
set(findobj('markerfacecolor','g'),'markerfacecolor','k');
figure(1);
set(get(gcf,'children'),'fontsize',12);

end;

%% Extra bonus scenario derived by thresholding LR04 stack. This generates a good fit
% if you allow the erosion rate and cover thickness to float as well. 

if 1;
close all;
ff = plotGISPdata(s)

mask = [1:7];

% range of threshold values
ds = linspace(3.2,3.4,20);

fits = ones(size(ds));
figure; semilogy(ds,fits,'ko','markerfacecolor','k','tag','ds'); drawnow;

% Scan across possible range. You can get a pretty good fit with 3.32-ish
% threshold, 10-ish m/Myr erosion, 150-ish g/cm2 cover. 
for a = 1:length(ds);
    X = [ds(a) 10 150];
    fits(a) = fit_LR04(X,s,p,m,mask,ff);
    set(findobj('tag','ds'),'ydata',fits); drawnow;
end;

end;

%% global optimization of LR04 stack threshold model

% Unfortunately this is terribly badly behaved because of the step-function behaviour
% of the threshold term. There are multiple minima and a downhill optimization
% doesn't work very well. It's actually easier to optimize this manually by
% spinning knobs -- threshold controls ratio, then change erosion rate to
% get absolute concentrations, then balance erosion rate with cover thickness to match
% depth profile. So you could maybe have the optimizer work in ratio vs.
% concentration space and it would work better.

if 1;

close all;
ff = plotGISPdata(s);

[optX,fmin] = fmincon(@(X) fit_LR04(X,s,p,m,mask,ff),[3.32 10 120],[],[],[],[],[2.8 0 0],[4 100 500],[],opts2)

end;




