Asian Option Arithmetic Turbunn Wakwman Form

function CallAsianOptionArithmeticTurbunnWakemanForm(S,K,r,sigma,q,T,t)
{
 sigmaC=sigma*sigma;
 M1=(Math.exp((r-q)*T)-Math.exp((r-q)*t))/((r-q)*(T-t));
 
 M2p1=2*Math.exp((2*(r-q)+sigmaC)*T)/((r-q+sigmaC)*(2*r-2*q+sigmaC)*(T-t)*(T-t));
 M2p2=2*Math.exp((2*(r-q)+sigmaC)*t)/((r-q)*(T-t)*(T-t))*((1/(2*(r-q)+sigmaC))-(Math.exp(((r-q))*T)/(r-q+sigmaC)));
 M2=M2p1+M2p2;
 b=ln(M1)/T;
 sigmaA=Math.sqrt((ln(M2)/T)-2*b);
 
 d1=(ln(S/K)+(b+0.5*sigmaC)*(T-t))/(sigmaA*Math.sqrt(T-t));
 d2=d1-(sigmaA*Math.sqrt(T-t));
 return S*Math.exp((b-r)*(T-t))*NCDF2(d1)-K*Math.exp(-r*(T-t))*NCDF2(d2);
}

function PutAsianOptionArithmeticTurbunnWakemanForm(S,K,r,sigma,q,T,t)
{
 sigmaC=sigma*sigma;
 M1=(Math.exp((r-q)*T)-Math.exp((r-q)*t))/((r-q)*(T-t));
 
 M2p1=2*Math.exp((2*(r-q)+sigmaC)*T)/((r-q+sigmaC)*(2*r-2*q+sigmaC)*(T-t)*(T-t));
 M2p2=(2*Math.exp((2*(r-q)+sigmaC)*t)/((r-q)*(T-t)*(T-t)))*((1/(2*(r-q)+sigmaC))-(Math.exp(((r-q))*T)/(r-q+sigmaC)));
 M2=M2p1+M2p2;
 b=ln(M1)/T;
 sigmaA=Math.sqrt((ln(M2)/T)-2*b);
 
 d1=(ln(S/K)+(b+0.5*sigmaC)*(T-t))/(sigmaA*Math.sqrt(T-t));
 d2=d1-(sigmaA*Math.sqrt(T-t));
 return K*Math.exp(-r*(T-t))*NCDF2(-d2)-S*Math.exp((b-r)*(T-t))*NCDF2(-d1);
}