function [K, Q] = KandQ(t, omega, c, d)
% KandQ -- compute the cov matrix and its inverse for the Fisher information
%
% Usage
% [K, Q] = KandQ(t, omega, c, d)
%
% Inputs
% t, omega, c, d estimated parameters
%
% Outputs
% K covariance matrix of H
% Q inverse of the covariance matrix
% Copyright (C) -- see DiscreteTFDs/Copyright
K = [...
1/16*(3+24*d^2*omega^2+24*d^4*c^2+16*d^4*omega^4+48*d^8*c^4+96*d^6*omega^2*c^2)/(d^4), -1/4+d^2*omega^2+3*d^4*c^2, -9/16*d^2+3/4*d^4*omega^2+15/4*d^6*c^2, 1/8*(-1+4*d^2*omega^2+12*d^4*c^2)/(d^4), -3/4*c-3*d^4*c^3-3*c*d^2*omega^2, 3/8*omega-1/2*d^2*omega^3-9/2*omega*d^4*c^2, 3*c*omega/d, 3*d^4*c*omega, -omega/d, -5/2*d*c ; ...
0, 3*d^4, 15/4*d^6, 3/2, -3*d^4*c, -3/2*d^4*omega, 0, 0, 0, 0 ; ...
0, 0, 105/16*d^8, 3/8*d^2, -15/4*d^6*c, -15/8*d^6*omega, 0, 0, 0, 0 ; ...
0, 0, 0, 9/4*1/(d^4), -3/2*c, -3/4*omega, 0, 0, 0, 0 ; ...
0, 0, 0, 0, d^2*omega^2+3/4+3*d^4*c^2, 3*d^4*c*omega, -omega/d, -3/2*d^4*omega, 0, d ; ...
0, 0, 0, 0, 0, 7/16*d^2+3/4*d^4*omega^2+15/4*d^6*c^2, -1/2*d*c, -15/4*d^6*c, 1/2*d, 0 ; ...
0, 0, 0, 0, 0, 0, 1/8*(4*d^2*omega^2+5+20*d^4*c^2)/(d^4), 0, -5/2*c, -5/4*omega ; ...
0, 0, 0, 0, 0, 0, 0, 15/4*d^6, 0, 0 ; ...
0, 0, 0, 0, 0, 0, 0, 0, 5/2, 0 ; ...
0, 0, 0, 0, 0, 0, 0, 0, 0, 39/8*d^2 ];
K = K + triu(K,1)';
Q = [...
1/24*(448*d^4*omega^4-1280*d^2*omega^2+2125)/(c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/288*(17920*d^6*omega^6+2700*d^4*c^2+16320*d^6*omega^2*c^2+27625-45376*d^4*omega^4+68360*d^2*omega^2-5376*d^8*c^2*omega^4)/(d^4*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/36*(-2112*d^4*omega^4+2970*d^2*omega^2+2125+896*d^6*omega^6)/(d^6*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/36*(-2112*d^4*omega^4+2970*d^2*omega^2+2125+896*d^6*omega^6)/(c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/6*(56*d^2*omega^2-95)/(c*(4*d^2*omega^2-5)), ...
-1/18*omega*(448*d^4*omega^4-1280*d^2*omega^2+2125)/(d^2*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
5/3*d^3*omega*(-59+8*d^2*omega^2)/(c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/9*omega*(224*d^4*omega^4-628*d^2*omega^2+1415)/(d^2*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/36*omega*(448*d^4*omega^4-3540*d^4*c^2-1280*d^2*omega^2+480*d^6*omega^2*c^2+2125)/(d*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/6*(8*d^2*omega^2-35)/(d*c*(4*d^2*omega^2-5)) ; ...
...
0, ...
1/3456*(-69120*d^10*c^4*omega^2+2144512*d^4*omega^4+162000*d^8*c^4+819000*d^4*c^2+1993680*d^2*omega^2-1582080*d^6*omega^6+1542656*d^8*c^2*omega^4+716800*d^8*omega^8+64512*d^12*c^4*omega^4-389120*d^10*c^2*omega^6+359125-493760*d^6*omega^2*c^2)/(d^8*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/432*(27625+91344*d^4*omega^4+53100*d^4*c^2+123610*d^2*omega^2-72832*d^6*omega^6+59520*d^8*c^2*omega^4+35840*d^8*omega^8-10752*d^10*c^2*omega^6-58920*d^6*omega^2*c^2)/(d^10*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/432*(27625+91344*d^4*omega^4+56700*d^4*c^2+123610*d^2*omega^2-72832*d^6*omega^6+61824*d^8*c^2*omega^4+35840*d^8*omega^8-10752*d^10*c^2*omega^6-64680*d^6*omega^2*c^2)/(d^4*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/72*(2240*d^4*omega^4-672*d^6*omega^2*c^2-3072*d^2*omega^2+540*d^4*c^2-1235)/((4*d^2*omega^2-5)*d^4*c), ...
1/216*omega*(27625-45376*d^4*omega^4+10700*d^4*c^2+68360*d^2*omega^2+17920*d^6*omega^6-256*d^8*c^2*omega^4+3520*d^6*omega^2*c^2)/(d^6*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-5/36*omega*(320*d^4*omega^4-468*d^4*c^2-767-2256*d^2*omega^2-96*d^6*omega^2*c^2)/(d*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/108*omega*(18395-22208*d^4*omega^4+7780*d^4*c^2+48436*d^2*omega^2+8960*d^6*omega^6-128*d^8*c^2*omega^4+3344*d^6*omega^2*c^2)/(d^6*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/432*omega*(-5760*d^10*c^4*omega^2+27625-45376*d^4*omega^4-28080*d^8*c^4-40120*d^4*c^2+68360*d^2*omega^2+17920*d^6*omega^6+15872*d^8*c^2*omega^4-124160*d^6*omega^2*c^2)/(d^5*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/72*(320*d^4*omega^4-96*d^6*omega^2*c^2-1296*d^2*omega^2-180*d^4*c^2-455)/((4*d^2*omega^2-5)*d^5*c) ; ...
...
0, ...
0, ...
1/54*(3828*d^4*omega^4+3600*d^4*c^2+2125+7220*d^2*omega^2-3328*d^6*omega^6+2304*d^8*c^2*omega^4+1792*d^8*omega^8-5760*d^6*omega^2*c^2)/(d^12*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/54*(3828*d^4*omega^4+3600*d^4*c^2+2125+7220*d^2*omega^2-3328*d^6*omega^6+2304*d^8*c^2*omega^4+1792*d^8*omega^8-5760*d^6*omega^2*c^2)/(d^6*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/9*(1+2*d^2*omega^2)*(56*d^2*omega^2-95)/((4*d^2*omega^2-5)*d^6*c), ...
-1/27*omega*(448*d^4*omega^4-1280*d^2*omega^2+2125)*(1+2*d^2*omega^2)/(d^8*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
10/9*omega*(-59+8*d^2*omega^2)*(1+2*d^2*omega^2)/(d^3*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-2/27*omega*(224*d^4*omega^4-628*d^2*omega^2+1415)*(1+2*d^2*omega^2)/(d^8*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/54*omega*(448*d^4*omega^4-3540*d^4*c^2-1280*d^2*omega^2+480*d^6*omega^2*c^2+2125)*(1+2*d^2*omega^2)/(d^7*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/9*(1+2*d^2*omega^2)*(8*d^2*omega^2-35)/((4*d^2*omega^2-5)*d^7*c) ; ...
...
0, ...
0, ...
0, ...
1/54*(3828*d^4*omega^4+4500*d^4*c^2+2125+7220*d^2*omega^2-3328*d^6*omega^6+2880*d^8*c^2*omega^4+1792*d^8*omega^8-7200*d^6*omega^2*c^2)/(c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/9*(1+2*d^2*omega^2)*(56*d^2*omega^2-95)/((4*d^2*omega^2-5)*c), ...
-1/27*omega*(448*d^4*omega^4-1280*d^2*omega^2+2125)*(1+2*d^2*omega^2)/(d^2*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
10/9*d^3*omega*(-59+8*d^2*omega^2)*(1+2*d^2*omega^2)/(c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-2/27*omega*(224*d^4*omega^4-628*d^2*omega^2+1415)*(1+2*d^2*omega^2)/(d^2*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/54*omega*(448*d^4*omega^4-3540*d^4*c^2-1280*d^2*omega^2+480*d^6*omega^2*c^2+2125)*(1+2*d^2*omega^2)/(d*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/9*(1+2*d^2*omega^2)*(8*d^2*omega^2-35)/((4*d^2*omega^2-5)*d*c) ; ...
...
0, ...
0, ...
0, ...
0, ...
14/3 , ...
-2/9*(56*d^2*omega^2-95)*omega/(d^2*c*(4*d^2*omega^2-5)), ...
20/3*d^3*omega/(4*d^2*omega^2-5), ...
-4/9*(28*d^2*omega^2-55)*omega/((4*d^2*omega^2-5)*d^2), ...
1/9*(56*d^2*omega^2+60*d^4*c^2-95)*omega/((4*d^2*omega^2-5)*d*c), ...
2/3*1/d ; ...
...
0, ...
0, ...
0, ...
0, ...
0, ...
2/27*(-1280*d^2*omega^4+1000*d^2*c^2+2125*omega^2+448*d^4*omega^6+640*d^6*c^2*omega^4-1600*d^4*omega^2*c^2)/(d^4*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-20/9*d*omega^2*(-59+8*d^2*omega^2)/(c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
4/27*(-628*d^2*omega^4+500*d^2*c^2+1415*omega^2+224*d^4*omega^6+320*d^6*c^2*omega^4-800*d^4*omega^2*c^2)/(d^4*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-1/27*(-1280*d^2*omega^4+400*d^2*c^2+2125*omega^2+448*d^4*omega^6+736*d^6*c^2*omega^4-4180*d^4*omega^2*c^2)/(d^3*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-2/9*(8*d^2*omega^2-35)*omega/((4*d^2*omega^2-5)*d^3*c) ; ...
...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
8/3*d^4*(15+37*d^2*omega^2)/(25+16*d^4*omega^4-40*d^2*omega^2), ...
-8/9*d*omega^2*(-221+20*d^2*omega^2)/(25+16*d^4*omega^4-40*d^2*omega^2), ...
2/9*d^2*(40*d^2*omega^4+180*d^2*c^2-295*omega^2+444*d^4*omega^2*c^2)/(c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
20/3*d^2*omega/(4*d^2*omega^2-5) ; ...
...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
4/27*(-800*d^6*omega^2*c^2+500*d^4*c^2+45-544*d^4*omega^4+320*d^8*c^2*omega^4+224*d^6*omega^6+1862*d^2*omega^2)/(d^6*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-2/27*(-628*d^2*omega^4+200*d^2*c^2+1415*omega^2+224*d^4*omega^6+368*d^6*c^2*omega^4-2972*d^4*omega^2*c^2)/(d^3*c*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
-4/9*(4*d^2*omega^2-25)*omega/((4*d^2*omega^2-5)*d^3) ; ...
...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
1/54*(5328*d^8*c^4*omega^2-1280*d^2*omega^4+2160*d^6*c^4+700*d^2*c^2+2125*omega^2+448*d^4*omega^6+1408*d^6*c^2*omega^4-8200*d^4*omega^2*c^2)/(d^2*c^2*(25+16*d^4*omega^4-40*d^2*omega^2)), ...
1/9*(8*d^2*omega^2+60*d^4*c^2-35)*omega/((4*d^2*omega^2-5)*c*d^2) ; ...
...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
0, ...
2/3*1/(d^2) ];
Q = Q + triu(Q,1)';