%enter the following
cin = 'momin';      %control file
out = 'momout';     %output file


%do not write anything below
NG=100;
LG=1000;
ntaper=[4 8 15];
fin=fopen(cin);
if fin==-1; disp('failure to open file'); disp(cin); quit; end;
fout=fopen(out,'w');
if fout==-1; disp('failure to open file'); disp(out); quit; end;
nburn=fscanf(fin,'%d',1);
nsave=fscanf(fin,'%d',1);
if nsave;
   apmfile=fscanf(fin,'%s',1); apmout=fopen(apmfile,'w');
   if apmout==-1; disp('failure to open file'); disp(apmfile); quit; end;
end;
hfile=fscanf(fin,'%s',1);
if nburn<0; nburn=0; end;

in=fopen(hfile);
if in==-1; fprintf(fout, 'failure to open file %s', hfile); quit; end;
nrun=0;
jf=fscanf(fin,'%d',1);
while (length(jf)~=0);
  err=0;
  if jf <= 0; 
    fprintf(fout, 'function of interest negative; skipped\n'); 
    err=1;
  else;
     frewind(in);
     niter=fscanf(in,'%d',1);   errcheck(niter,fout,0);
     nf=fscanf(in,'%d',1);      errcheck(nf,fout,0);
     if (jf > nf);
        fprintf(fout,'function of interest %3d greater than ',jf);
        fprintf(fout,'limit of %3d , skipped\n',nf);
        err=1;
     elseif (nburn>niter);
        fprintf(fout,'%3d burn in requested but only ',nburn);
        fprintf(fout,'%3d avaliable', niter);
        quit;
     elseif (nf>LG);
        fprintf(fout,'file contains %3d functions per record, ',nf);
        fprintf(fout,'exceeds storage capacity of %3d ', LG);
        quit;
     end;
  end; %(jf<=0)
 
  if (~err);
     nrun=nrun+1;
     for iburn=1:nburn;
        iter=fscanf(in,'%d',1); errcheck(iter,fout,iter);
        wtlog=fscanf(in,'%lf',1); errcheck(wtlog,fout,iter);
        dummy1=fscanf(in,'%lf',1); errcheck(dummy1,fout,iter);
        dummy2=fscanf(in,'%lf',1); errcheck(dummy2,fout,iter);
        for i=1:nf;
           g(i)=fscanf(in,'%lf',1); errcheck(g(i),fout,iter);
        end;
%       fprintf('%20.8e \n',dummy1); 
     end;
     nuse=niter-nburn;
     if (nuse<NG);
        fprintf(fout, 'minimum %3d analysis iterations are required; ',NG);
        fprintf(fout, '%3d avaliable \n', nuse);
        term;
     end;
     ns=floor(nuse/NG);
     nuse=ns*NG;
     if (nrun==1); header1(fout,niter,nburn,nuse); end;

% form sufficiency statistics for all the subsequent work

     td=0; tn=0; tdd=0; tnn=0; tdn=0; tvar=0;
     for ig=1:NG;
        gd=0; gn=0; gdd=0; gdn=0; gnn=0; gvar=0;
        for is=1:ns;
           iter=fscanf(in,'%d',1); errcheck(iter,fout,iter);
           wtlog=fscanf(in,'%lf',1); errcheck(wtlog,fout,iter);
           dummy1=fscanf(in,'%lf',1); errcheck(dummy1,fout,iter);
           dummy2=fscanf(in,'%lf',1); errcheck(dummy2,fout,iter);
           for i=1:nf;
              g(i)=fscanf(in,'%lf',1); errcheck(g(i),fout,iter);
           end;
           ad=0;
           if (wtlog > -700);
              ad=exp(wtlog);
           end;
           an=ad*g(jf);
           gd=gd+ad;
           gn=gn+an;
           gdn=gdn+ad*an;
           gdd=gdd+ad*ad;
           gnn=gnn+an*an;
           gvar=gvar + an*g(jf);
        end; % for is
        td=td+gd;
        tn=tn+gn;
        tdn=tdn+gdn;
        tdd=tdd+gdd;
        tnn=tnn+gnn;
        tvar=tvar+gvar;

        cn(ig)=gn/ns;
        cd(ig)=gd/ns;
        cdn(ig)=gdn/ns; 
        cdd(ig)=gdd/ns;
        cnn(ig)=gnn/ns;
        cvar(ig)=gvar/ns;
     end; %for ig
     eg=tn/td;
     varg=tvar/td -eg^2;
     sdg=-1; if (varg>0); sdg=sqrt(varg); end;
     fprintf(fout, '\nFunction%3d Posterior mean:%16.7e', jf, eg);
     fprintf(fout, ' Posterior st dev:%16.7e\n\n',sdg); 
     if nsave;
        fprintf(apmout, '%4d%6d%16.7e%16.7e\n', jf, nuse, eg, sdg); 
     end;

%numerical standard error assuming no serial correlation

     varnum=(tnn-2*eg*tdn+tdd*eg^2)/(td^2);
     sdnum=-1; if (varnum>0); sdnum=sqrt(varnum); end; 
     header2(fout);
     fprintf(fout,'%20.7e %15.6f\n', sdnum, varg/(nuse*varnum));
     if nsave;
        fprintf(apmout, '%16.7e  ', sdnum); 
     end;

%get autocovariance of grouped means
     barn=tn/nuse;
     bard=td/nuse;
     for ig=1:NG;
       cn(ig)=cn(ig)-barn;
       cd(ig)=cd(ig)-bard;
     end;
     for lag=0:NG-1;
        ann=0; add=0; and=0; adn=0;
        for ig=lag+1:NG;
           ann=ann+cn(ig)*cn(ig-lag);
           add=add+cd(ig)*cd(ig-lag);
           and=and+cn(ig)*cd(ig-lag);
           adn=adn+cd(ig)*cd(ig-lag);
        end; %ig
 % index 0 not allowed, lag+1 stands for lag
        rnn(lag+1)=ann/NG;
        rdd(lag+1)=add/NG;
        rnd(lag+1)=and/NG;
        rdn(lag+1)=adn/NG;
     end; %lag

% numerical standard error with tapered autocovariance functions
     for mm=1:3;
        m=ntaper(mm);
        am=m;
        snn=rnn(1); sdd=rdd(1); snd=rnd(1);
        for lag=1:m-1;
           att=1-lag/am;
           snn=snn+2*att*rnn(lag+1);
           sdd=sdd+2*att*rdd(lag+1);
           snd=snd+att*(rnd(lag+1) + rnd(lag+1)); 
        end; %lag
        varnum=ns*nuse*(snn-2*eg*snd+sdd*eg^2)/(td^2);
        sdnum=-1; if (varnum>0); sdnum=sqrt(varnum); end;
        c='%';
        fprintf(fout, 'Autocovariance function tapered to %3d %1s', m,c);
        fprintf(fout, ' %20.7e %15.6f\n', sdnum, varg/(nuse*varnum)); 
        if nsave;
           fprintf(apmout, '%16.7e  ', sdnum); 
        end;
      end; %mm
      if nsave;
         fprintf(apmout, '\n', sdnum); 
      end;
  end; % (~err)
  jf=fscanf(fin,'%d',1);    %read the next function number
end;  %length(jf)


%if the program exits the while loop without reaching EOF it must
%be the case that there is I/O error with the input for jf 
if ~feof(fin); 
  fprintf(fout, 'error on input for function number, program stops'); 
  quit;
end;
if nrun==0;
  fprintf(fout, 'error on input for function number, program stops'); 
  quit;
end;
fclose(fin);
fclose(fout);
quit;