function [llf]=mmlemods(p,pid,p2snam,perts,gg)
%                         USED BY MMLE.M
 % CREATES MODELS FOR COMPILED FORTRAN GHES.EXE ROUTINE TO CALC THE
 % GRADIENT AND HESSIAN OF THE COST FUNCTION. MODELS STORED IN TOGHES.MAT.
 % GHES.EXE READS TOGHES.MAT FOR THE SYSTEM MATRICES, UYDAT.MAT FOR THE
 % TIME HISTORIES, AND WRITES FROMGHES.MAT WHICH CONTAINS GRAD AND HES.

%define global variable

global uydata
global inovt
global wersum
global grad
global rr
global e0
global grade
global rowinq
global diag_rr


dnorm=0;% USED TO CHECK IF D MATRICES ALL ZERO, TO SAVE COMPUTATIONS IN GHES

p2snamp=[p2snam,'(p)']; %                             ACCEPT ANY NAME FOR P2SS
[a,phi,gam,c,d,q,x0,dt,rowinq]=eval(p2snamp);c0=c; %        GET NOMINAL SYSTEM
npid=length(pid);[nm,nx]=size(c);[ndp,xx]=size(uydata); %       GET DIMENSIONS
[nx,nu]=size(gam);num=nu+nm;nnorm=ndp;
grade=ones(nx,npid); %                                      DIMENSION MATRICES
proc=sum(diag(q))>0; %  USER SETS q TO ZERO FOR NO KALMAN GAIN
 %------------------YEST,INOV,RR,LLF FOR NOMINAL P----------------------------
if proc==0,
      k=zeros(nx,nm);
else
      k=lqe(a,eye(nx),c,q*q'/dt/dt,gg)*dt; % EQN 30,MAINE & ILIFF, pp588-579
end%                                   SIAM J. Appl. Math. Vol41, No3, Dec81
 %
phinom=phi*(eye(nx)-k*c);gamnom=[gam-phi*k*d,phi*k];
     dnom=[-d,eye(nm)];
dnorm=dnorm+norm(d,1);
inovt=dlsim(phinom,gamnom,-c,dnom,uydata,x0);
 %
 %
rr=inovt'*inovt/nnorm; %            INNOVATIONS AND SAMPLE COVARIANCE
rtggi=chol(inv(gg)); %    OUTPUT VECTOR WEIGHTING FOR COST FUNCTION
inovt=inovt*rtggi';%yest=yest*rtggi'; %             WEIGHT  YEST,INOVT
 %-------------------------------------------------------------------
e0=diag(k*c); %                    FOR CONSTRAINING DIAG(E0) < EYE
wersum=sum(sum(inovt.*inovt))/nnorm/nm;%        WEIGHTED ERROR SUM
llf=ndp/2*(wersum*nm-2*log(prod(diag(rtggi))));%
 %----------------------------NOW GET GRAD, HES, OF COST WRT PARAMS
psav=p;
s=[phinom,gamnom;rtggi*[c,d],zeros(nm)];xpid=x0';  % STACK NOMINAL PLANT

k=zeros(nx,nm);
for jj=1:npid; %        COMPUTE AND STACK KF MATRICES FOR ALL PERTURBATIONS
  p(pid(jj))=p(pid(jj))+perts(pid(jj));
  [a,phi,gam,c,d,q,x0,dt,rowinq]=eval(p2snamp); %      GET PERTURBED SYSTEM
  if proc ~ 0; %                      KALMAN FILTER IF PROCESS NOISE
 %      k=lqe(a,eye(nx),c,q*q'/dt/dt,gg)*dt; %   CONT LQE FOR DISC! (SEE REF)
    [zzv,zzd]=eig([a' c'/gg*c;q*q'/dt/dt -a]);% THESE 5 LINES EQUIV TO LQE
    zzd=diag(zzd);[zzd,index]=sort(real(zzd));
    chi=zzv(1:nx,index(1:nx));
    lambda=zzv(nx+1:2*nx,index(1:nx));zzs=-real(lambda/chi);
    k=gg\c*zzs;k=k'*dt;
  end;
  dnorm=dnorm+norm(d,1); %
  s=[s;[phi*(eye(nx)-k*c),[gam-phi*k*d,phi*k];rtggi*[c,d],zeros(nm)]];
  grade(:,jj)=(diag(k*c)-e0)/perts(pid(jj)); %        GRADIENT OF DIAG(E=K*C)
 %
  xpid=[xpid x0']; %                STACK INITIAL CONDITIONS IN ROWS
  p=psav;
  fprintf(' ')
end;
fprintf('.')
 %------------------------------------------interface to ghes.exe
s=s';
size([phi gamnom]);abcol=ans(2);cdcol=nx+nu*(dnorm>0);
rtggit=rtggi';
toghes=[npid,nx,nm,abcol,cdcol,num,-7,-8,-9,-10,perts(pid),xpid(:)',...
rtggit(:)',s(:)', zeros(1,50)];
save toghes toghes
%----------------------------------------- end

