function [T] = random_3stages(d, z_fix,z_pc)

% create 3 stage ice history (simpler than Balco et al.) Stage 1 (S1) = post LGM
% thin ice cover, stage 2 = exposure/thinner than present, stage 3 = re-thickening at present
% thickness.

% generate random start point 10 ka - 6 ka
Tmax = randi([6000, 10000]);
rand_Ts = sort(randi([60, Tmax],1,2));% Now choose two random times between Tmax and 60 a (TMA shows ice config the same in 1962)
T.all_Ts = [0 rand_Ts(1) rand_Ts rand_Ts(2) Tmax]; % concatenate (6 nodes)

% do initial ice thickness
Z1 = randi([50, 200]); % random initial ice thickness between 200 and 50 m (i.e. some iterations allow small amount of muon production)

if z_fix == 1 % not allowing ice thickness to be totally free but fixing % of runs to be 0 thickness
            % S2 thickness is between 0 and present ice thickness but want % of runs to be at 0 m
            m = 100 / z_pc;
        r = randi(m); % random 1,2,3,4 etc to define as %
        if r == 1
        Z2 = 0; % fix at 0 for ~set % of runs 
        else
         Z2 = rand(1) .* d.ice_thickness; % others its jsut random
        end
else
    Z2 = rand(1) .* d.ice_thickness; % all iterations are random
end
    
Z3 = d.ice_thickness; % present ice thickness
T.all_Zs = [Z3 Z3 Z2 Z2 Z1 Z1]; % concatenate thickness 6 nodes

T.ice_surf_eleTmax_T1 = round(d.surf_elv + Z1); % get ice surf elevation for 1st stage for PR calcs (not used)
