function E = meade_bsa_obj_synth(IN)
%
%Meade triangular dislocation objective function 
%For use with Gibbs sampling algorithm 
%Data are digitized seismic horizons from the BSA
%
%USAGE E = meade_bsa_obj(IN);
%
%--------------------------------
%SETUP - SPECIFY FIXED PARAMETERS
%--------------------------------
l_fix=1e6;%length
dy_fix=l_fix/2;%upper left hand corner
st_fix=180;%strike
ss_fix=0;%strike-slip component
%
%----------------------------------
%SETUP - CREATE TRIANGLES FOR MEADE
%----------------------------------
%function below takes input dislocation geometry and creates
%properly oriented trianguler dislocation elements
%tdng_jw=make_tdng_jw(dx,dy,l,w,d,strike,dip)
%tdng_jw=make_tdng_jw(IN.m(1),IN.m(2),IN.m(3),IN.m(4),IN.m(5),IN.m(6),IN.m(7));
tdng_jw=make_tdng_jw(IN.m(1),dy_fix,l_fix,IN.m(4),IN.m(5),st_fix,IN.m(7));
%
%triangle vertex ordering
vpos_tri=tdng_jw.vpos;vpos_tri(3,:)=-vpos_tri(3,:);
tv1=[tdng_jw.tv(1,1) tdng_jw.tv(2,1) tdng_jw.tv(3,1)];
tv2=[tdng_jw.tv(1,2) tdng_jw.tv(2,2) tdng_jw.tv(3,2)];
ts=0;%tensile slip displacement set to zero
%
%----------
%RUN MEADE
%----------
%CalcTriDisps(sx, sy, sz, x, y, z, pr, ss, ts, ds)
%where sx,sy,sz are observation points
%and x,y,z are triangle vertex locations
%and ss,ts,ds are imposed displacements
x1=vpos_tri(1,tv1);
y1=vpos_tri(2,tv1);
z1=vpos_tri(3,tv1);
U1=CalcTriDisps(IN.x',IN.y',IN.z', x1, y1, z1, IN.nu, ss_fix, ts, IN.m(9));
x2=vpos_tri(1,tv2);
y2=vpos_tri(2,tv2);
z2=vpos_tri(3,tv2);
U2=CalcTriDisps(IN.x',IN.y',IN.z', x2, y2, z2, IN.nu, ss_fix, ts, IN.m(9));
%superpose displacement fields for triangular elements
U3z=U1.z+U2.z;
U3x=U1.x+U2.x;
%Z model points output from Meade code above (kilometers)
%Model=[IN.x' -IN.z'+U3z];
%X/Z model points output from Meade code above (kilometers)
%spline to regional x position
%hard-wired number of layers
%n=length(U3z);
%s1=spline(IN.x(1:n/2)'+U3x(1:n/2),-IN.z(1:n/2)'+U3z(1:n/2),IN.x(1:n/2)');
%s2=spline(IN.x(n/2+1:end)'+U3x(n/2+1:end),-IN.z(n/2+1:end)'+U3z(n/2+1:end),IN.x(n/2+1:end)');
%Model=[IN.x' [s1;s2]];
%loop over the number of layers (works with two...test with more)
n=IN.length_layers{1};
s=spline(IN.x(1:n)'+U3x(1:n),-IN.z(1:n)'+U3z(1:n),IN.x(1:n)');
if IN.num_layers>1;
    for i=2:IN.num_layers
        str=length(s)+1;
        stp=str+IN.length_layers{i}-1;
        ss=spline(IN.x(str:stp)'+U3x(str:stp),-IN.z(str:stp)'+U3z(str:stp),IN.x(str:stp)');
        s=[s;ss];%#ok<AGROW>
    end
end
Model=[IN.x' s];
%-------
%ERRORS
%-------
resid=IN.Obs-Model; %full 2-d vector
resid=reshape(resid',2*length(IN.Obs),1);
E=resid'*IN.Cinv*resid./(length(Model)-5);%(d-Gm)'*C-1*(d-Gm)
%no covariance matrix for no-noise synthetic case
%E=resid'*resid./(length(Model)-length(IN.m));%(d-Gm)'*(d-Gm)
%E=resid'*resid./(length(Model)-IN.free_params);%(d-Gm)'*(d-Gm)
%L2=sqrt(resid'*IN.OBJ.Cinv*resid)./length(Model);