function [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,...
   nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,...
   valueuv,matrix_bar,criterion)
%OBJECTIVEFUNCTION_C - compute the misfit of a model
%
%   this routine updates the value of the spatial distribution of viscosity, calculates the
%   new velocity field and compares it to the observed velocity field (J)
%
%   Usage:
%      [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion)

B_new=B+search_scalar*direction;
%@@@ B_bar=matrix_bar*B_new;

B_bar=(B_new(index(:,1))+B_new(index(:,2))+B_new(index(:,3)))/3;
%### nu_bar=10^14*ones(nel,1);

%###AK
nu_bar=viscosity(index,nel,alpha,beta,u_objective,v_objective,B_bar,glen_coeff);
%###

converged_yet=0;
convergence_count=1;
loop=0;


while (~converged_yet)
   clear u v
   
   
   nu_bar_uu=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
   nu_bar_vv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
   nu_bar_uv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
   
   if (loop>100)
		u=u_old; v=v_old;
      break;
   end;
   nu_bar_oldvalue=nu_bar;
   loop=loop+1;
   
   Duu=sparse(rowDshort,colDshort,nu_bar_uu.*valueuu,nods,nods);
   Duu=Duu+triu(Duu,1)';
   
   Dvv=sparse(rowDshort,colDshort,nu_bar_vv.*valuevv,nods,nods);
   Dvv=Dvv+triu(Dvv,1)';
   
   F=[Duu sparse(rowD,colD,nu_bar_uv.*valueuv,nods,nods)
      sparse(colD,rowD,nu_bar_uv.*valueuv,nods,nods) Dvv]; 
   
   
   Rhs_parsed=P*(Rhs - F*specified_velocity);
   Fprime=P*F*P';
   
   
   solution=Fprime\Rhs_parsed;
   solution=P'*solution + specified_velocity;
   u=solution(1:nods);
   v=solution(nods+1:2*nods);
   nu_bar=viscosity(index,nel,alpha,beta, u,v,B_bar,glen_coeff);
   
   %Test for direct shooting convergence
   if convergence_count>1,

       dug=[u; v]-[u_old; v_old];
       ndug=norm(dug,2);
       ug=[u_old;v_old];
       nug=norm(ug,2);

       change=ndug/nug;

       %Figure out if viscosity converged
       if change>criterion,
           converged_yet=0;
       else
           converged_yet=1;
       end
   else
       converged_yet=0;
   end

   u_old=u;
   v_old=v;

convergence_count=convergence_count+1;
 
end 
%disp(['# of elements over criterion: ' num2str(sum(test))]);

%@@@AK save F_file F
save F_file F 

Du=[(u-vx_obs).*weighting
   (v-vy_obs).*weighting];


J=Du'*S*Du;
