clear all
close all

runid = 'chan_M1e4';
pmcrit = 0.3;

n = 6;
m = 76;
t = 51;


for i=1:11
    trimfile = ['trim-',runid,'_',num2str(i-1),'.dat'];
    TRIM      = vs_use(trimfile,'quiet');
    pm1(i)    = vs_get(TRIM,'map-sed-series',{t},'MUDFRAC',{n,m},'quiet');
    Esand1(i) = vs_get(TRIM,'map-sed-series',{t},'EFLUX',{n,m,1},'quiet');
    Emud1(i)  = vs_get(TRIM,'map-sed-series',{t},'EFLUX',{n,m,2},'quiet');     
    f1(i)     = pm1(i)/(1-pm1(i));
    if (1-pm1(i))>0
        Es1(i)    = Esand1(i)/(1-pm1(i))/1600;
    else 
        Es1(i)    = 0;
    end
    Es1(i)    = f1(i)*Es1(i);
    if pm(i)>0
        Em1(i)    = Emud1(i)/pm1(i)/500;
    else 
        Em1(i)    = 0;
    end  
    if pm1(i)<=pmcrit
        sandmud1(i) = Es1(i);
    else
        sandmud1(i) = Em1(i)*(Es1(i)/Em1(i))^((1-pm1(i))/(1-pmcrit));
    end
end

figure(10)
hold on
grid on
plot(pm1,Em1,'g')
plot(pm1,Es1,'k')
plot(pm1,sandmud1,'k--')
% ylim([0 5e-6])


for i=1:11
    trimfile = ['trim-',runid,'_pm03_',num2str(i-1),'.dat'];
    TRIM  = vs_use(trimfile,'quiet');
    pm2(i)     = vs_get(TRIM,'map-sed-series',{t},'MUDFRAC',{n,m},'quiet');
    Esand2(i) = vs_get(TRIM,'map-sed-series',{t},'EFLUX',{n,m,1},'quiet');
    Emud2(i)  = vs_get(TRIM,'map-sed-series',{t},'EFLUX',{n,m,2},'quiet');      
    f2(i)     = pm2(i)/(1-pm2(i));
    if pm2(i)>0
        Em2(i)    = Emud2(i)/pm2(i)/500;
    else 
        Em2(i)    = 0;
    end
    if (1-pm(i))>0
        Es2(i)    = Esand2(i)/(1-pm2(i))/1600;
    else 
        Es2(i)    = 0;
    end
    Es2(i)    = f2(i)*Es2(i);
    if pm2(i)<=pmcrit
        sandmud2(i) = Es2(i);
    else
        sandmud2(i) = Em1(i)*(Es2(i)/Em1(i))^((1-pm2(i))/(1-pmcrit));
    end
end

plot(pm2,Em2,'r','Linewidth',2)
plot(pm2,Es2,'b')
% plot(pm,sandmud2,'b--')