program dkkmake(input,output);
{ compute columns for DK }
const
  hdim=29;
  edim=59;
  lmax=3;
  vdim=200;
  rho_val='6.0E+0';
  ncolumns=29;

%include 'types.h';

var
  normalize: boolean;
  i,j: integer;
  beta2,rho,cov,sigma,tau,st: scalar;
  u,umax,lnormal: ltuple;
  vkappa,ve,v1,v2,vt1,vt2: vector;
  vv: vtuple;
  t: text;

%include 'reps.i';
%include 'scalar.i';
%include 'vector.i';
%include 'linear.i';
%include 'kkmap.i';

begin { dkkmake }
  sinit;
  ssset(rho_val,rho);
  sread('beta2.s',beta2);
  sread('cov.s',cov);
  vread('approx.v',rho,vt1);
  vread('kkappa.v',rho,vkappa);
  vread('kkball.v',rho,v1);
  writeln;

  kkinit(vt1,u);         { approx.v to be checked by dkkcheck }
  writeln;

  writeln('dkkmake: initializing tangent map DK');
  normalize:=true;
  ustar(v1,u);
  dkkinit(u,v1,ve,vv);
  vsprod(u[0],vv[0],vt1);
  vscale(beta2,vt1,v2);         { v2 = unpacked and scaled v1 }
  writeln;

  writeln('dkkmake: computing columns of DK');
  rewrite(t,'dkkt.m');
  writeln(t,ncolumns);

  for i:=0 to ncolumns do begin
    writeln('column',i:3);
    vzero(vt1);
    if i<ncolumns then vt1[i]:=sone                 { order i }
    else lvho(ncolumns,rho,sone,vt1);    { higher order terms }

    dkk(i,u,vv,ve,v2,vt1,vt2);

    for j:=0 to ncolumns-1 do begin
      writeln(t,rst(vt2[j].l));
      writeln(t,rst(vt2[j].u));
    end;
    vlho(ncolumns,vt2,rho,st);
    writeln(t,rst(st.l));
    writeln(t,rst(st.u));
    writeln;
  end;
  flush(t);

  sdone('dkkmake');
end { dkkmake }.