function [pomega]=p(omega,t,q)

cq=pi/(q-1)*(gamma((1/(q-1)-0.5))/gamma(1/(q-1)))^2;
bta=cq^((1-q)/(3-q))*(t.*(2-q)*(3-q)).^(-2/(3-q));
zt=(t.*cq*(2-q)*(3-q))^(1/(3-q));

pomega=1./zt.*(1+omega.^2.*bta.*(q-1)).^(-1/(q-1));