function [Call, Put] = COSmethod(Strike, params, r, Spot, T, model, JT)
% Computation of the cumulants of the log-asset price
[c1, c2, c4] = GetCumulants(r, params, T, JT);

% Truncation range
L = 12;
a = c1-L*sqrt(c2+sqrt(c4));
b = c1+L*sqrt(c2+sqrt(c4));

% Pricing european options using Fourier-cosine series
Npow = 10;
N = 2^Npow;     % Number of elements in the series
k = 0:N-1;

x = log(Spot./Strike);
m = length(x);
ExpMatrix = zeros(m,N);
for i = 1:m
    ExpMatrix(i,:) = exp(1i*pi*k*(x(i,1)-a)/(b-a));
end

w = ones(1,N);      % weights
w(1) = 0.5;     % First term of the series is weighted by 1/2

% Computation of Put prices
c = a;
d = 0;
add1 = cos(k*pi*(d-a)/(b-a))*exp(d) - cos(k*pi*(c-a)/(b-a))*exp(c); 
add2 = ((k*pi)/(b-a)) .* (sin(k*pi*(d-a)/(b-a))*exp(d) - sin(k*pi*(c-a)/(b-a))*exp(c));
fatt = 1./(1+((k*pi)/(b-a)).^2);
chi  = fatt .* (add1+add2);

psi = (sin(k*pi*(d-a)/(b-a)) - sin(k*pi*(c-a)/(b-a))) .* ((b-a)./(k*pi));
psi(1) = d-c;

U_k = 2/(b-a)*(-chi+psi);

u = k*pi/(b-a);
phi = CharFun(params, r, T, u, model, JT);

Put = zeros(m,1);
for i = 1:m
    Put(i) = exp(-r*T)*Strike(i)*real(w.*phi.*ExpMatrix(i,:))*U_k';
end

% Call prices are obtained through Put-Call parity
Call = Put+Spot-Strike.*exp(-r*T);

end