% Global profile found by: sd=0.03+0.23*(tanh(2:-5.5/48:-3.5)+1)

sd=[ ...
    0.5000 0.5000
    0.4817 0.5000
    0.4796 0.5000
    0.4771 0.5000
    0.4738 0.5000
    0.4699 0.5000
    0.4649 0.5000
    0.4589 0.5000
    0.4516 0.5000
    0.4427 0.5000
    0.4321 0.5000
    0.4194 0.5000
    0.4046 0.5000
    0.3876 0.5000
    0.3682 0.5000
    0.3466 0.5000
    0.3230 0.5000
    0.2980 0.5000
    0.2720 0.5000
    0.2456 0.5000
    0.2197 0.5000
    0.1948 0.5000
    0.1714 0.5000
    0.1500 0.5000
    0.1308 0.5000
    0.1139 0.5000
    0.0993 0.5000
    0.0869 0.5000
    0.0764 0.5000
    0.0677 0.5000
    0.0605 0.5000
    0.0546 0.5000
    0.0498 0.5000
    0.0458 0.5000
    0.0427 0.5000
    0.0401 0.5000
    0.0381 0.5000
    0.0365 0.5000
    0.0352 0.5000
    0.0341 0.5000
    0.0333 0.5000
    0.0326 0.5000
    0.0321 0.5000
    0.0317 0.5000
    0.0313 0.5000
    0.0310 0.5000
    0.0308 0.5000
    0.0307 0.5000
    0.0305 0.5000
    0.0304 0.5000
    0.1000 0.5000
    0.2200 0.5000
    ];

%sigbot=[37.800+0.47 39.0000];    % cguo: 0.47 comes from eosben07_rho(2000,3.5,36)-eosben07_rho(2000,3.5,35.4)
sigbot=[37.800+0.78 39.0000];    % cguo: 0.78 comes from eosben07_rho(2000,3.5,36.4)-eosben07_rho(2000,3.5,35.4)

i_index=[    1  142];
j_index=[   27  282];

nzm=size(sd,1)+1;
ndm=length(i_index);
ndt=length(sigbot);
if ndm>ndt
  sigbot(ndt+1:ndm)=sigbot(ndt);
  sd(:,ndt+1:ndm)=sd(:,ndt)*ones(1,ndm-ndt);
end

sig_n=zeros(nzm,ndm);
sig_n(end,1:ndm)=sigbot(1:ndm);

for k=nzm-1:-1:1
  sig_n(k,1:ndm)=sig_n(k+1,1:ndm)-sd(k,1:ndm);
end
[(1:nzm)' sig_n]

fid=fopen('sigref.txt','w');
fprintf(fid,'%7d%7d\n',nzm,ndm);
fprintf(fid,'%7d',i_index);
fprintf(fid,'\n');
fprintf(fid,'%7d',j_index);
fprintf(fid,'\n');
for k=1:nzm
  fprintf(fid,'%7.3f',sig_n(k,:));
  fprintf(fid,'\n');
end
fclose(fid);
