clf
%17 parameters = vp(1:3), vs(1:3), rho(1:3), qp(1:3), qs(1:3), thickness(1:2)
    [obs_data, dt, rayp, npts] = ReadSac('stack_ALI.sac');
    obs_data = obs_data ./max(obs_data);
    thickness = [5 30 0];
    vp = [5.0 6.5 8.0];
    rho = [2.8 2.8 3.3];
    vs = (1/sqrt(3)).*vp;    
    qs = 100.*vs;
    qp = 200.*vs;

    if (rayp < 0), rayp = 0.05; end
    theta = ((180/pi)*asin(rayp*vp(2)));

    model = [thickness; vp; vs; rho; qp; qs];
    [time, RF] = RecFun(theta, model, dt, npts);

    plot(time, RF,'b--')
    hold on
    N = max(size(time));
    
    % line up model and observed maximum value
    [xmax, ymax1] = max(obs_data);
    [xmax, ymax2] = max(RF);
    diff = ymax1 - ymax2;
    
    % plot observed receiver function (red) and calculate L1 error
    plot(time, obs_data(1+diff:N+diff),'r');
    errst = sum(abs(obs_data(1+diff:N+diff) - (RF')));
    err0=errst;
    
       thicknessst(1)=thickness(1);
       thicknessst(2)=thickness(2);
       thicknessst(3)=0;
       thickness0(1)=thickness(1);
       thickness0(2)=thickness(2);
       thickness0(3)=0;
for it=1:300

it

thickness(1)=2.5+5*rand;
thickness(2)=15.+30.*rand;
vp(1)=2.5+5.*rand;
vp(2)=3.25+6.5*rand;
vp(3)=7.8+0.4*rand;
vs = (1/sqrt(3)).*vp;    
    qs = 100.*vs;
    qp = 200.*vs;
%

    if (rayp < 0), rayp = 0.05; end
    theta = ((180/pi)*asin(rayp*vp(2)));

    model = [thickness; vp; vs; rho; qp; qs];
    [time, RF] = RecFun(theta, model, dt, npts);

    N = max(size(time));
   
    % line up model and observed maximum value
    [xmax, ymax1] = max(obs_data);
    [xmax, ymax2] = max(RF);
    diff = ymax1 - ymax2;
   
    % plot observed receiver function (red) and calculate L1 error
%   plot(time, obs_data(1+diff:N+diff),'r');
    err(it) = sum(abs(obs_data(1+diff:N+diff) - (RF')));
    if (err(it)<err0), 
       err0=err(it);
       it0=it;
       vp0=vp;
       vs0=vs;
       rho0=rho;
       qs0=qs;
       qp0=qp;
    end;

    % print out crustal parameters
    fprintf ('------------------------------\n');
    fprintf('Moho depth: %6.0f km    Error: %4.1f\n', thickness(1) + thickness(2), err(it));
    if (vp(1) > vp(2)) 
        fprintf('****Warning: low velocity zone!\n');
        fprintf('****Upper crust Vp > lower crust Vp\n');
    end
end;

    if (rayp < 0), rayp = 0.05; end
    theta = ((180/pi)*asin(rayp*vp(2)));

    model = [thickness; vp0; vs0; rho0; qp0; qs0];
    [time, RF] = RecFun(theta, model, dt, npts);

    N = max(size(time));
   
    plot(time, RF,'g')

        fprintf ('Starting model (SM) err: %4.0f\n', errst);
    for i = 1:3
        fprintf ('BM Layer %d depth: %4.1f Vp: %6.2f Vs: %6.2f', i, thickness(i), vp(i), vs(i));
        fprintf ('BM rho: %5.2f Qp: %4.0f Qs: %4.0f\n', rho(i), qp(i),qs(i));
    end

        fprintf ('Best Model (BM) model err: %4.0f\n', err(it0));
    for i = 1:3
        fprintf ('BM Layer %d depth: %4.1f Vp: %6.2f Vs: %6.2f', i, thickness(i), vp0(i), vs0(i));
        fprintf ('BM rho: %5.2f Qp: %4.0f Qs: %4.0f\n', rho0(i), qp0(i),qs0(i));
    end