out ccon;
%off nat;
on rounded;

procedure delta(ii,jj);
begin 
if ii eq jj then return 1 else return 0; 
end;

n:=15;m:=6;
operator a,k,r,rm,f;

for all ii such that ii<0 or ii>m-1 let k(ii)=0;
for all ii such that ii<0 or ii>n-1 let a(ii)=0;
for all ii such that ii<0 or ii>n-1 let r(ii)=0;
for all ii such that ii<0 or ii>n-1 let rm(ii)=0;

for ii:=-10:20 do write(k(ii));
for ii:=-10:20 do write(a(ii));
for ii:=-10:20 do write(r(ii));
for ii:=-10:20 do write(rm(ii));

%-------------A-----------
for ii:=0:(n-1) do a(ii):=10; a(7):=1000; a(9):=1000;
%-------------centered psf kernel-----------
for ii:=0:(m-1) do k(ii):=0; 
k(0):=0.2; k(1):=0.3; k(2):=0.5; 
k(3):=0.5; k(4):=0.3; k(5):=0.2; 

%-------------convolution procedure-----------
%psf kernel dim = m
%------------------
procedure cconv();
begin

for tt:=0:(n-1) do <<

if tt > (m - 1) then write( r(tt) := for ii:=0:(m-1) sum a(tt-ii)*k(ii) )
		else write(r(tt):=0);
>>;
end;
%---------------------------------------------
%-------------convolution procedure-----------
%psf kernel dim = n
%---------------------------------------------

procedure cconvn();
begin

for tt:=0:(n-1) do <<

if tt > (m - 1) then write( r(tt) := for ii:=0:(n-1) sum a(tt-ii)*k(ii) )
		else write(r(tt):=0);
>>;
end;
%---------------------------------------------
procedure cconvm();
begin

for tt:=0:(n-1) do <<

%if tt > (m - 1) then 

write( rm(tt) := for pp:=0:(n-1) sum 
(for ii:=0:(n-1) sum k(ii)*delta(pp,tt-ii))*A(pp)

     )


		%else write(r(tt):=0);
>>;
end;
%---------------------------------------------
for ii:=0:(m-1) do a(ii):=0;
%---------------------------------------------

cconv();
cconvn();

for ii:=0:(m-1) do write("k(",ii,"):=",k(ii) );
for ii:=0:(n-1) do write("a(",ii,"):=",a(ii) );
for ii:=0:(n-1) do write("r(",ii,"):=",r(ii) );

cconvm();
for ii:=0:(n-1) do write("rm(",ii,"):=",r(ii) );

%for ii:=0:(n-1) do write(rr(ii):=r(ii+m/2) );



for ii:=0:(n-1) do write( ii, "   ", r(ii) - rm(ii) );



matrix AA(1,n);
for ii:=1:n do AA(1,ii):=a(ii-1);
AA;

matrix RR(1,n);
for ii:=1:n do rr(1,ii):=r(ii-1);
RR;

matrix RRM(1,n);
for ii:=1:n do rrm(1,ii):=rm(ii-1);
RRM;


matrix KK(n,n);
for tt:=1:n do for pp:=1:n do
kk(tt,pp) := for ii:=0:(n-1) sum k(ii)*delta(pp-1,tt-1-ii);

KK;

tp(rr) - kk*tp(aa);

rr(1,8):=205;
ikk:=kk**-1;
ikk*tp(rr) - tp(aa);




shut ccon;
bye;
end;