// gives two decimal places
function twod(x) {z=Math.round(100*x)/100; return z}

// function to find a root of.  This function is, approximately speaking, supply - demand in the intermediate market (raised to a convenient power)
function f(ty) {temp1=0.0;temp2=0.0;
for (ig=1; ig<=n; ig=ig+1) {temp1=temp1+(Math.pow(psi(ty,stemp[ig],rtemp[ig]),beta))*k[ig];temp2=temp2+(Math.pow(chi(ty,stemp[ig],rtemp[ig]),eta))*g[ig]};
vf= Math.pow(temp1,1+a*eta)-Math.pow(temp2,1+a*beta);return vf}

// Define psi, chi functions 
// psi and chi are not the same functions from the appendix, have 1-t-psi, t-xi instead
//This prevents a literal zero, which gives conniptions when raised to a power, but doesn't affect actual zeros based on market shares.
function psi(th,st,rt) {
vpsi=Math.max(0.00000001,1-th-(((1-th)/beta)*((th*(st-rt)/eta)+(a*st*(1-rt)))/((a*(1-st)*(1-rt))+((1-th)*(1-rt)/beta)+(th*(1-st)/eta)))); return vpsi}

function chi(th,st,rt) {
vchi=Math.max(0.00000001,th-((th/eta)*(((1-th)*(rt-st)/beta)+(a*rt*(1-st)))/((a*(1-st)*(1-rt))+((1-th)*(1-rt)/beta)+(th*(1-st)/eta)))); return vchi}

// this function reads the inputs
function cal() {
alpha=document.mhi.alpha.value;alpha=alpha*1;
beta=document.mhi.beta.value;beta=beta*1;
eta=document.mhi.eta.value;eta=eta*1;
n=document.mhi.n.value;n=n*1;
theta=document.mhi.theta.value;theta=theta*1;
type=document.mhi.type.value;type=type*1;
sb[1]=document.mhi.sb1.value; rb[1]=document.mhi.rb1.value;sb[1]=sb[1]*1;rb[1]=rb[1]*1;
sb[2]=document.mhi.sb2.value; rb[2]=document.mhi.rb2.value;sb[2]=sb[2]*1;rb[2]=rb[2]*1;
sb[3]=document.mhi.sb3.value; rb[3]=document.mhi.rb3.value;sb[3]=sb[3]*1;rb[3]=rb[3]*1;
sb[4]=document.mhi.sb4.value; rb[4]=document.mhi.rb4.value;sb[4]=sb[4]*1;rb[4]=rb[4]*1;
sb[5]=document.mhi.sb5.value; rb[5]=document.mhi.rb5.value;sb[5]=sb[5]*1;rb[5]=rb[5]*1;
sb[6]=document.mhi.sb6.value; rb[6]=document.mhi.rb6.value;sb[6]=sb[6]*1;rb[6]=rb[6]*1;
sb[7]=document.mhi.sb7.value; rb[7]=document.mhi.rb7.value;sb[7]=sb[7]*1;rb[7]=rb[7]*1;
sb[8]=document.mhi.sb8.value; rb[8]=document.mhi.rb8.value;sb[8]=sb[8]*1;rb[8]=rb[8]*1;
sb[9]=document.mhi.sb9.value; rb[9]=document.mhi.rb9.value;sb[9]=sb[9]*1;rb[9]=rb[9]*1;
sb[10]=document.mhi.sb10.value; rb[10]=document.mhi.rb10.value;sb[10]=sb[10]*1;rb[10]=rb[10]*1;
sb[11]=document.mhi.sb11.value; rb[11]=document.mhi.rb11.value;sb[11]=sb[11]*1;rb[11]=rb[11]*1;
sb[12]=document.mhi.sb12.value; rb[12]=document.mhi.rb12.value;sb[12]=sb[12]*1;rb[12]=rb[12]*1;
sb[13]=document.mhi.sb13.value; rb[13]=document.mhi.rb13.value;sb[13]=sb[13]*1;rb[13]=rb[13]*1;
sb[14]=document.mhi.sb14.value; rb[14]=document.mhi.rb14.value;sb[14]=sb[14]*1;rb[14]=rb[14]*1;
sb[15]=document.mhi.sb15.value; rb[15]=document.mhi.rb15.value;sb[15]=sb[15]*1;rb[15]=rb[15]*1;
sb[16]=document.mhi.sb16.value; rb[16]=document.mhi.rb16.value;sb[16]=sb[16]*1;rb[16]=rb[16]*1;
sb[17]=document.mhi.sb17.value; rb[17]=document.mhi.rb17.value;sb[17]=sb[17]*1;rb[17]=rb[17]*1;
sb[18]=document.mhi.sb18.value; rb[18]=document.mhi.rb18.value;sb[18]=sb[18]*1;rb[18]=rb[18]*1;
sb[19]=document.mhi.sb19.value; rb[19]=document.mhi.rb19.value;sb[19]=sb[19]*1;rb[19]=rb[19]*1;
sb[20]=document.mhi.sb20.value; rb[20]=document.mhi.rb20.value;sb[20]=sb[20]*1;rb[20]=rb[20]*1;
}

// this function rewrites the input variables and is called if one of the inputs in inadmissible
function displayfix() {
document.mhi.alpha.value=alpha;
document.mhi.beta.value=beta;
document.mhi.eta.value=eta;
document.mhi.n.value=n;
document.mhi.theta.value=theta;
document.mhi.type.value=type;
document.mhi.sb1.value=sb[1]; document.mhi.rb1.value=rb[1];
document.mhi.sb2.value=sb[2]; document.mhi.rb2.value=rb[2];
document.mhi.sb3.value=sb[3]; document.mhi.rb3.value=rb[3];
document.mhi.sb4.value=sb[4]; document.mhi.rb4.value=rb[4];
document.mhi.sb5.value=sb[5]; document.mhi.rb5.value=rb[5];
document.mhi.sb6.value=sb[6]; document.mhi.rb6.value=rb[6];
document.mhi.sb7.value=sb[7]; document.mhi.rb7.value=rb[7];
document.mhi.sb8.value=sb[8]; document.mhi.rb8.value=rb[8];
document.mhi.sb9.value=sb[9]; document.mhi.rb9.value=rb[9];
document.mhi.sb10.value=sb[10]; document.mhi.rb10.value=rb[10];
}

// Note that sb, rb are in percent, and capital values are in absolute levels, while s, r are not in percent

// the main computational program
function comput() {

errmessage="Messages:\n";stp=0.5;dfix=0; //errmessage is the messages file, stp is step size, and dfix=1 means some input was illegal
document.mhi.SOCA.style.backgroundColor="#FFFFFF";document.mhi.SOCA.value="";
document.mhi.SOCB.style.backgroundColor="#FFFFFF";document.mhi.SOCB.value="";

// Error checking
if(type!=1&&type!=2&&type!=3) {type=1;dfix=1;errmessage=errmessage+"\nError: Illegal Merger Type, reset to 1"}
if(type==1) {errmessage=errmessage+"\nThe mergers combines both upstream and downstream assets of firms 1 and 2"};
if(type==2) {errmessage=errmessage+"\nThe merger combines downstream assets of firms 1 and 2 but spins off 2's upstream assets."};
if(type==3) {errmessage=errmessage+"\nThe merger combines upstream assets of firms 1 and 2 but spins off 2's downstream assets."};
if(!(n>=3)){errmessage=errmessage+"\nLow, negative or nonnumeric n reset to 3";n=3;dfix=1;}
n=Math.round(n);
if(!(alpha>=0||alpha<=0)){alpha=3;dfix=1;errmessage=errmessage+"\nError: Nonnumeric alpha encountered"}
if(alpha<=0) {a=0} else {a=1/alpha}
if(alpha<=1&&alpha>0) {errmessage=errmessage+"\nWarning: With inelastic demand, monopoly solution is q=0, and solution may not exist."}
if(alpha<=0) {errmessage=errmessage+"\nDemand elasticity set to infinity.  This produces the special case"}
if(!(beta>0)) {beta=1;dfix=1;errmessage=errmessage+"\nError: Downstream cost elasticities must be positive numbers."}
if(!(eta>0)) {eta=1;dfix=1;errmessage=errmessage+"\nError: Upstream cost elasticities must be positive numbers."}
if(!(theta>0&&theta<1)){theta=0.5;dfix=1;errmessage=errmessage+"\nError: By definition, 0<theta<1."}

// calculate input share totals and put negative input shares to zero
sbsum=0.0; rbsum=0.0;
for (i=1; i<=n; i=i+1) {if(!(sb[i]>=0)) {dfix=1;sb[i]=0}; if(!(rb[i]>=0)) {dfix=1;rb[i]=0};sbsum=sbsum+sb[i];rbsum=rbsum+rb[i]}

//check problems with inputs
if(sbsum<=0) {dfix=1;errmessage=errmessage+"\nError: Negative or zero downstream market shares.";sb[1]=20;sb[2]=30;sb[3]=50;sbsum=100;n=Math.max(n,3)}
if(rbsum<=0) {dfix=1;errmessage=errmessage+"\nError: Negative or zero upstream market shares.";rb[1]=20;rb[2]=30;rb[3]=50;rbsum=100;n=Math.max(n,3)}
if(sbsum!=100) {dfix=1;for (i=1; i<=n; i=i+1) {sb[i]=twod(100*sb[i]/sbsum)};errmessage=errmessage+"\nError: Downstream shares don't sum to 100% and were rescaled"}
if(rbsum!=100) {dfix=1;for (i=1; i<=n; i=i+1) {rb[i]=twod(100*rb[i]/rbsum)};errmessage=errmessage+"\nError: Upstream shares don't sum to 100% and were rescaled"}

//zero out shares for inputs beyond n
overnsum=0;
for (i=n+1; i<=10; i=i+1) {overnsum=overnsum+sb[i];overnsum=overnsum+sb[i];sb[i]=0;rb[i]=0;s[i]=0;r[i]=0;kb[i]=0;gb[i]=0;k[i]=0;g[i]=0;stemp[i]=0;rtemp[i]=0;};
if(overnsum>0) {dfix=1};

if(!(dfix==0)) {displayfix();errmessage=errmessage+"\nRewriting input values to eliminate errors and inconsistencies"};

//Check second order conditions
socb="Passed";
errmessage=errmessage+"\n \n...Checking Second Order Conditions...";
for (i=1; i<=n; i=i+1)
{term1=d2k(sb[i]/100,rb[i]/100,a,(1-theta)/beta,theta/eta,theta);
term2=d2g(sb[i]/100,rb[i]/100,a,(1-theta)/beta,theta/eta,theta);
term3=crossterm(sb[i]/100,rb[i]/100,a,(1-theta)/beta,theta/eta,theta);

if (term1<=0&&term2<=0&&term3>=0) {errmessage=errmessage+"\nSOCs satisfied for firm ";errmessage=errmessage+i}
else {socb="Failed";errmessage=errmessage+"\n\nSOCs failed for firm ";errmessage=errmessage+i;
errmessage=errmessage+"\nThis set of inputs is not consistent with profit maximization!!";
errmessage=errmessage+"\nterm 1 = ";errmessage=errmessage+term1;
errmessage=errmessage+"\nterm 2 = ";errmessage=errmessage+term2;
errmessage=errmessage+"\nterm 3 = ";errmessage=errmessage+term3}

} //End SOC checking
document.mhi.SOCB.value=socb;
document.mhi.error.value=errmessage;
if (socb=="Failed") {document.mhi.SOCB.style.backgroundColor="#FF0000"} else {document.mhi.SOCB.style.backgroundColor="#00FF00"}

// Begin Computations

for (i=1; i<=n; i=i+1) {s[i]=sb[i]/100;r[i]=rb[i]/100;   }

ch=1;for (i=1; i<=n; i=i+1) {ch=Math.min(ch,psi(theta,s[i],r[i]));ch=Math.min(ch,chi(theta,s[i],r[i]))}
if(!(ch>0)) {errmessage=errmessage+"\nError: Likely failure of existence, try increasing alpha"}

document.mhi.error.value=errmessage;

// Compute capital levels k, g for q=1, k is the same notation as the paper, g is gamma
// then store starting values
// stemp, rtemp will be share estimates to iterate to a solution of the new shares  
ksum=0; gsum=0;
for (i=1; i<=n; i=i+1) {kb[i]=(Math.pow(psi(theta,s[i],r[i]),-beta))*s[i];gb[i]=(Math.pow(chi(theta,s[i],r[i]),-eta))*r[i];k[i]=kb[i];g[i]=gb[i];stemp[i]=s[i];rtemp[i]=r[i];ksum=ksum+k[i];gsum=gsum+g[i]}

// Perform Merger on capitals, estimate shares 
if(type==1) {k[1]=k[1]+k[2];k[2]=0;g[1]=g[1]+g[2];g[2]=0;s[1]=s[1]+s[2];s[2]=0;r[1]=r[1]+r[2];r[2]=0;stemp[1]=stemp[1]+stemp[2];stemp[2]=0;rtemp[1]=rtemp[1]+rtemp[2];rtemp[2]=0}
if(type==2) {k[1]=k[1]+k[2];k[2]=0;s[1]=s[1]+s[2];s[2]=0;stemp[1]=stemp[1]+stemp[2];stemp[2]=0}
if(type==3) {g[1]=g[1]+g[2];g[2]=0;r[1]=r[1]+r[2];r[2]=0;rtemp[1]=rtemp[1]+rtemp[2];rtemp[2]=0}

// end input

// This do loop computes shares, iterate to find a solution 
for (iter=1; iter<=2000; iter=iter+1) {
  
xleft=0;xright=1;
// Find f[x]==0

if(f(xleft)*f(xright)>0) {xstar=theta; errmessage=errmessage+"\nApparently no solution for f==0.  Phooey."}

if(f(xleft)<=0&&f(xright)>=0) {
for (im=1; im<=30; im=im+1) {xtest=(xleft+xright)/2; if(f(xtest)>0) {xright=xtest} else {xleft=xtest} }
xstar=(xleft+xright)/2}

if(f(xleft)>=0&&f(xright)<=0) {
for (im=1; im<=30; im=im+1) {xtest=(xleft+xright)/2; if(f(xtest)>0) {xleft=xtest} else {xright=xtest} }
xstar=(xleft+xright)/2}

thetatemp=xstar;
  
// Now compute shares that should prevail given the new values of q and theta 

sqerror=0;
tempq=0.0;for (iq=1; iq<=n; iq=iq+1) {tempq=tempq+(Math.pow(psi(thetatemp,stemp[iq],rtemp[iq]),beta)*k[iq])}
qtemp=Math.pow(tempq,1/(1+a*beta));

ssum=0;rsum=0;
for (iw=1; iw<=n; iw=iw+1) {
s[iw]=((1-stp)*s[iw])+(stp*Math.pow(qtemp,-a*beta-1)*Math.pow(psi(thetatemp,stemp[iw],rtemp[iw]),beta)*k[iw]);
r[iw]=((1-stp)*r[iw])+(stp*Math.pow(qtemp,-a*eta-1)*Math.pow(chi(thetatemp,stemp[iw],rtemp[iw]),eta)*g[iw]);
ssum=ssum+s[iw];rsum=rsum+r[iw]}
for (iw=1; iw<=n; iw=iw+1) {
s[iw]=s[iw]/ssum;r[iw]=r[iw]/rsum;
sqerror=sqerror+((stemp[iw]-s[iw])*(stemp[iw]-s[iw]))+((rtemp[iw]-r[iw])*(rtemp[iw]-r[iw]));
stemp[iw]=s[iw];rtemp[iw]=r[iw]  }

// if we have a solution, the calculation should give the same shares we started with 
bb=iter;bsq=sqerror;if(Math.sqrt(sqerror)<error) {bb=iter;break}
if(iter%25==0){stp=stp*2/3;errmessage=errmessage+"\nslow convergence, reducing step size"}

 } // End of Iter loop

q=qtemp;
t=thetatemp;

// The efficient quantity and intermediate good price 

ql=0.0;qr=Math.max(1,ksum);
for(iq=1; iq<100; iq=iq+1) {
qq=(ql+qr)/2;
if(((Math.pow(qq,a+(1/beta)))/(Math.pow(ksum,1/beta)))+((Math.pow(qq,a+(1/eta)))/(Math.pow(gsum,1/eta)))<1) {ql=qq} else {qr=qq}  }

qe=qq;

pe=Math.pow(qe,a+(1/eta))/Math.pow(gsum,1/eta);

MHIold=0.0;MHI=0.0;
for(iu=1; iu<=n; iu=iu+1) {stemp[iu]=sb[iu]/100;rtemp[iu]=rb[iu]/100;
MHIold=MHIold+((theta*(1-theta)*((stemp[iu]-rtemp[iu])*(stemp[iu]-rtemp[iu]))/(eta*beta))+(a*(1-theta)*(1-rtemp[iu])*(stemp[iu]*stemp[iu])/beta)+(a*theta*(1-stemp[iu])*(rtemp[iu]*rtemp[iu])/eta))/((a*(1-stemp[iu])*(1-rtemp[iu]))+((1-theta)*(1-rtemp[iu])/beta)+(theta*(1-stemp[iu])/eta));
MHI= MHI+((t*(1-t)*((s[iu]-r[iu])*(s[iu]-r[iu]))/(eta*beta))+(a*(1-t)*(1-r[iu])*(s[iu]*s[iu])/beta)+(a*t*(1-s[iu])*(r[iu]*r[iu])/eta))/((a*(1-s[iu])*(1-r[iu]))+((1-t)*(1-r[iu])/beta)+(t*(1-s[iu])/eta))
}

// End computations, start output 

if(!(Math.sqrt(sqerror)<error)) {
errmessage=errmessage+"\nError: Failure to converge.  Iterations= ";errmessage=errmessage+bb;errmessage=errmessage+" Estimated error is ";errmessage=errmessage+sqerror}
else {errmessage=errmessage+"\nConvergence achieved at iteration ";errmessage=errmessage+bb}

for(i=1; i<=n; i=i+1) {k[i]=k[i]/ksum;g[i]=g[i]/gsum;kb[i]=kb[i]/ksum;gb[i]=gb[i]/gsum}

document.mhi.MHIold.style.backgroundColor="#FFFFFF";
document.mhi.qold.style.backgroundColor="#FFFFFF";
document.mhi.MHI.style.backgroundColor="#FFFFFF";
document.mhi.q.style.backgroundColor="#FFFFFF"
if(MHIold>=1) {errmessage=errmessage+"\n\nLikely nonsense answer.  Initial MHI>=100%.  Probably non-existence of solution.";document.mhi.MHIold.style.backgroundColor="#FF3333"}
if(1<0.0001*qe) {errmessage=errmessage+"\n\nLikely nonsense answer.  Initial q very near zero.  Probably non-existence of solution.";document.mhi.qold.style.backgroundColor="#FF3333"}
if(MHI>=1) {errmessage=errmessage+"\n\nLikely nonsense answer.  MHI>=100%.  Probably non-existence of solution.";document.mhi.MHI.style.backgroundColor="#FF3333"}
if(q<0.0001*qe) {errmessage=errmessage+"\n\nLikely nonsense answer.  Calculated q very near zero.  Probably non-existence of solution.";document.mhi.q.style.backgroundColor="#FF3333"}

// write output

document.mhi.error.value=errmessage;
document.mhi.kb1.value=twod(100*kb[1]);document.mhi.gb1.value=twod(100*gb[1]);
document.mhi.kb2.value=twod(100*kb[2]);document.mhi.gb2.value=twod(100*gb[2]);
document.mhi.kb3.value=twod(100*kb[3]);document.mhi.gb3.value=twod(100*gb[3]);
document.mhi.kb4.value=twod(100*kb[4]);document.mhi.gb4.value=twod(100*gb[4]);
document.mhi.kb5.value=twod(100*kb[5]);document.mhi.gb5.value=twod(100*gb[5]);
document.mhi.kb6.value=twod(100*kb[6]);document.mhi.gb6.value=twod(100*gb[6]);
document.mhi.kb7.value=twod(100*kb[7]);document.mhi.gb7.value=twod(100*gb[7]);
document.mhi.kb8.value=twod(100*kb[8]);document.mhi.gb8.value=twod(100*gb[8]);
document.mhi.kb9.value=twod(100*kb[9]);document.mhi.gb9.value=twod(100*gb[9]);
document.mhi.kb10.value=twod(100*kb[10]);document.mhi.gb10.value=twod(100*gb[10]);
document.mhi.kb11.value=twod(100*kb[11]);document.mhi.gb11.value=twod(100*gb[11]);
document.mhi.kb12.value=twod(100*kb[12]);document.mhi.gb12.value=twod(100*gb[12]);
document.mhi.kb13.value=twod(100*kb[13]);document.mhi.gb13.value=twod(100*gb[13]);
document.mhi.kb14.value=twod(100*kb[14]);document.mhi.gb14.value=twod(100*gb[14]);
document.mhi.kb15.value=twod(100*kb[15]);document.mhi.gb15.value=twod(100*gb[15]);
document.mhi.kb16.value=twod(100*kb[16]);document.mhi.gb16.value=twod(100*gb[16]);
document.mhi.kb17.value=twod(100*kb[17]);document.mhi.gb17.value=twod(100*gb[17]);
document.mhi.kb18.value=twod(100*kb[18]);document.mhi.gb18.value=twod(100*gb[18]);
document.mhi.kb19.value=twod(100*kb[19]);document.mhi.gb19.value=twod(100*gb[19]);
document.mhi.kb20.value=twod(100*kb[20]);document.mhi.gb20.value=twod(100*gb[20]);
document.mhi.error.value=errmessage;

document.mhi.s1.value=twod(100*s[1]);document.mhi.r1.value=twod(100*r[1]);document.mhi.k1.value=twod(100*k[1]);document.mhi.g1.value=twod(100*g[1]);
document.mhi.s2.value=twod(100*s[2]);document.mhi.r2.value=twod(100*r[2]);document.mhi.k2.value=twod(100*k[2]);document.mhi.g2.value=twod(100*g[2]);
document.mhi.s3.value=twod(100*s[3]);document.mhi.r3.value=twod(100*r[3]);document.mhi.k3.value=twod(100*k[3]);document.mhi.g3.value=twod(100*g[3]);
document.mhi.s4.value=twod(100*s[4]);document.mhi.r4.value=twod(100*r[4]);document.mhi.k4.value=twod(100*k[4]);document.mhi.g4.value=twod(100*g[4]);
document.mhi.s5.value=twod(100*s[5]);document.mhi.r5.value=twod(100*r[5]);document.mhi.k5.value=twod(100*k[5]);document.mhi.g5.value=twod(100*g[5]);
document.mhi.s6.value=twod(100*s[6]);document.mhi.r6.value=twod(100*r[6]);document.mhi.k6.value=twod(100*k[6]);document.mhi.g6.value=twod(100*g[6]);
document.mhi.s7.value=twod(100*s[7]);document.mhi.r7.value=twod(100*r[7]);document.mhi.k7.value=twod(100*k[7]);document.mhi.g7.value=twod(100*g[7]);
document.mhi.s8.value=twod(100*s[8]);document.mhi.r8.value=twod(100*r[8]);document.mhi.k8.value=twod(100*k[8]);document.mhi.g8.value=twod(100*g[8]);
document.mhi.s9.value=twod(100*s[9]);document.mhi.r9.value=twod(100*r[9]);document.mhi.k9.value=twod(100*k[9]);document.mhi.g9.value=twod(100*g[9]);
document.mhi.s10.value=twod(100*s[10]);document.mhi.r10.value=twod(100*r[10]);document.mhi.k10.value=twod(100*k[10]);document.mhi.g10.value=twod(100*g[10]);
document.mhi.s11.value=twod(100*s[11]);document.mhi.r11.value=twod(100*r[11]);document.mhi.k11.value=twod(100*k[11]);document.mhi.g11.value=twod(100*g[11]);
document.mhi.s12.value=twod(100*s[12]);document.mhi.r12.value=twod(100*r[12]);document.mhi.k12.value=twod(100*k[12]);document.mhi.g12.value=twod(100*g[12]);
document.mhi.s13.value=twod(100*s[13]);document.mhi.r13.value=twod(100*r[13]);document.mhi.k13.value=twod(100*k[13]);document.mhi.g13.value=twod(100*g[13]);
document.mhi.s14.value=twod(100*s[14]);document.mhi.r14.value=twod(100*r[14]);document.mhi.k14.value=twod(100*k[14]);document.mhi.g14.value=twod(100*g[14]);
document.mhi.s15.value=twod(100*s[15]);document.mhi.r15.value=twod(100*r[15]);document.mhi.k15.value=twod(100*k[15]);document.mhi.g15.value=twod(100*g[15]);
document.mhi.s16.value=twod(100*s[16]);document.mhi.r16.value=twod(100*r[16]);document.mhi.k16.value=twod(100*k[16]);document.mhi.g16.value=twod(100*g[16]);
document.mhi.s17.value=twod(100*s[17]);document.mhi.r17.value=twod(100*r[17]);document.mhi.k17.value=twod(100*k[17]);document.mhi.g17.value=twod(100*g[17]);
document.mhi.s18.value=twod(100*s[18]);document.mhi.r18.value=twod(100*r[18]);document.mhi.k18.value=twod(100*k[18]);document.mhi.g18.value=twod(100*g[18]);
document.mhi.s19.value=twod(100*s[19]);document.mhi.r19.value=twod(100*r[19]);document.mhi.k19.value=twod(100*k[19]);document.mhi.g19.value=twod(100*g[19]);
document.mhi.s20.value=twod(100*s[20]);document.mhi.r20.value=twod(100*r[20]);document.mhi.k20.value=twod(100*k[20]);document.mhi.g20.value=twod(100*g[20]);

document.mhi.displaytheta.value=twod(100*theta);
document.mhi.t.value=twod(100*t);
document.mhi.perctheta.value=twod(100.0*(t-theta)/theta);
document.mhi.pe.value=twod(100*pe);
document.mhi.qold.value=twod(100.0/qe);
document.mhi.q.value=twod(100.0*q/qe);
document.mhi.percq.value=twod(100.0*(q-1));
document.mhi.MHIold.value=twod(100.0*MHIold);
document.mhi.MHI.value=twod(100.0*MHI);

//Check second order conditions
soca="Passed";
errmessage=errmessage+"\n \n...Checking second order conditions for calculated solution...";
for (i=1; i<=n; i=i+1)
{term1=d2k(s[i],r[i],a,(1-t)/beta,t/eta,t);
term2=d2g(s[i],r[i],a,(1-t)/beta,t/eta,t);
term3=crossterm(s[i],r[i],a,(1-t)/beta,t/eta,t);

if (term1<=0&&term2<=0&&term3>=0) {errmessage=errmessage+"\nSOCs satisfied for firm ";errmessage=errmessage+i}
else {soca="Failed";errmessage=errmessage+"\n\nSOCs failed for firm ";errmessage=errmessage+i;
errmessage=errmessage+"\nterm 1 = ";errmessage=errmessage+term1;
errmessage=errmessage+"\nterm 2 = ";errmessage=errmessage+term2;
errmessage=errmessage+"\nterm 3 = ";errmessage=errmessage+term3}

} //End SOC checking
document.mhi.SOCA.value=soca;
document.mhi.error.value=errmessage;
if (soca=="Failed") {document.mhi.SOCA.style.backgroundColor="#FF0000"} else {document.mhi.SOCA.style.backgroundColor="#00FF00"}

}  //End of computing section

function d2k(si,ri,A,B,C,t) {
A2=A*A;A3=A2*A;A4=A3*A;B2=B*B;B3=B2*B;B4=B3*B;C2=C*C;C3=C2*C;C4=C3*C;
t2=t*t;t3=t2*t;t4=t3*t;si2=si*si;si3=si2*si;ri2=ri*ri;ri3=ri2*ri;
calc=A*B*C2*ri*si;
calc=calc + B2*C2*ri*si;
calc=calc - B2*C2*si2;
calc=calc - A*B*C2*ri*si2;
calc=calc - A3*t;
calc=calc - 3*A2*B*t;
calc=calc - 3*A*B2*t;
calc=calc - B3*t;
calc=calc - 3*A2*C*t;
calc=calc - 6*A*B*C*t;
calc=calc - 3*B2*C*t;
calc=calc - 3*A*C2*t;
calc=calc - 3*B*C2*t;
calc=calc - C3*t;
calc=calc + A3*ri*t;
calc=calc + 3*A2*B*ri*t;
calc=calc + 3*A*B2*ri*t;
calc=calc + B3*ri*t;
calc=calc + 2*A2*C*ri*t;
calc=calc + 4*A*B*C*ri*t;
calc=calc - A2*B*C*ri*t;
calc=calc + 2*B2*C*ri*t;
calc=calc - 2*A*B2*C*ri*t;
calc=calc - B3*C*ri*t;
calc=calc + A*C2*ri*t;
calc=calc + B*C2*ri*t;
calc=calc - 2*A*B*C2*ri*t;
calc=calc - 2*B2*C2*ri*t;
calc=calc - B*C3*ri*t;
calc=calc + A3*si*t;
calc=calc + 2*A2*B*si*t;
calc=calc + A3*B*si*t;
calc=calc + A*B2*si*t;
calc=calc + 2*A2*B2*si*t;
calc=calc + A*B3*si*t;
calc=calc + 3*A2*C*si*t;
calc=calc + 4*A*B*C*si*t;
calc=calc + 3*A2*B*C*si*t;
calc=calc + B2*C*si*t;
calc=calc + 4*A*B2*C*si*t;
calc=calc + B3*C*si*t;
calc=calc + 3*A*C2*si*t;
calc=calc + 2*B*C2*si*t;
calc=calc + 3*A*B*C2*si*t;
calc=calc + 2*B2*C2*si*t;
calc=calc + C3*si*t;
calc=calc + B*C3*si*t;
calc=calc - A3*ri*si*t;
calc=calc - 2*A2*B*ri*si*t;
calc=calc - A3*B*ri*si*t;
calc=calc - A*B2*ri*si*t;
calc=calc - 2*A2*B2*ri*si*t;
calc=calc - A*B3*ri*si*t;
calc=calc - A2*C*ri*si*t;
calc=calc - A*B*C*ri*si*t;
calc=calc + A*C2*ri*si*t;
calc=calc + B*C2*ri*si*t;
calc=calc + B2*C2*ri*si*t;
calc=calc + C3*ri*si*t;
calc=calc + B*C3*ri*si*t;
calc=calc + A2*B*si2*t;
calc=calc - A3*B*si2*t;
calc=calc + A*B2*si2*t;
calc=calc - A2*B2*si2*t;
calc=calc + 2*A*B*C*si2*t;
calc=calc - 3*A2*B*C*si2*t;
calc=calc + B2*C*si2*t;
calc=calc - 4*A*B2*C*si2*t;
calc=calc + B*C2*si2*t;
calc=calc - 3*A*B*C2*si2*t;
calc=calc - B2*C2*si2*t;
calc=calc - B*C3*si2*t;
calc=calc - A2*B*ri*si2*t;
calc=calc + A3*B*ri*si2*t;
calc=calc - A*B2*ri*si2*t;
calc=calc + A2*B2*ri*si2*t;
calc=calc - A2*C*ri*si2*t;
calc=calc - 3*A*B*C*ri*si2*t;
calc=calc + A2*B*C*ri*si2*t;
calc=calc - B2*C*ri*si2*t;
calc=calc + 2*A*B2*C*ri*si2*t;
calc=calc - 2*A*C2*ri*si2*t;
calc=calc - 2*B*C2*ri*si2*t;
calc=calc + 2*A*B*C2*ri*si2*t;
calc=calc - C3*ri*si2*t;
calc=calc + A3*t2;
calc=calc + 3*A2*B*t2;
calc=calc + 3*A*B2*t2;
calc=calc + B3*t2;
calc=calc + 3*A2*C*t2;
calc=calc + 6*A*B*C*t2;
calc=calc + 3*B2*C*t2;
calc=calc + 3*A*C2*t2;
calc=calc + 3*B*C2*t2;
calc=calc + C3*t2;
calc=calc - A3*ri*t2;
calc=calc - 3*A2*B*ri*t2;
calc=calc - 3*A*B2*ri*t2;
calc=calc - B3*ri*t2;
calc=calc - 2*A2*C*ri*t2;
calc=calc - 4*A*B*C*ri*t2;
calc=calc - 2*B2*C*ri*t2;
calc=calc - A*C2*ri*t2;
calc=calc - B*C2*ri*t2;
calc=calc - A3*si*t2;
calc=calc - 2*A2*B*si*t2;
calc=calc - A*B2*si*t2;
calc=calc - 3*A2*C*si*t2;
calc=calc - 4*A*B*C*si*t2;
calc=calc - B2*C*si*t2;
calc=calc - 3*A*C2*si*t2;
calc=calc - 2*B*C2*si*t2;
calc=calc - C3*si*t2;
calc=calc + A3*ri*si*t2;
calc=calc + 2*A2*B*ri*si*t2;
calc=calc + A*B2*ri*si*t2;
calc=calc + A2*C*ri*si*t2;
calc=calc + A*B*C*ri*si*t2;
calc=calc - A2*B*C*ri*si*t2;
calc=calc - A*C2*ri*si*t2;
calc=calc - B*C2*ri*si*t2;
calc=calc - C3*ri*si*t2;
calc=calc - A2*B*si2*t2;
calc=calc - A*B2*si2*t2;
calc=calc - A2*B2*si2*t2;
calc=calc - 2*A*B*C*si2*t2;
calc=calc - B2*C*si2*t2;
calc=calc - B*C2*si2*t2;
calc=calc + A2*B*ri*si2*t2;
calc=calc + A*B2*ri*si2*t2;
calc=calc + A2*B2*ri*si2*t2;
calc=calc + A2*C*ri*si2*t2;
calc=calc + 3*A*B*C*ri*si2*t2;
calc=calc + A2*B*C*ri*si2*t2;
calc=calc + B2*C*ri*si2*t2;
calc=calc + 2*A*C2*ri*si2*t2;
calc=calc + 2*B*C2*ri*si2*t2;
calc=calc + C3*ri*si2*t2; return calc}

function d2g(si,ri,A,B,C,t) {
A2=A*A;A3=A2*A;A4=A3*A;B2=B*B;B3=B2*B;B4=B3*B;C2=C*C;C3=C2*C;C4=C3*C;
t2=t*t;t3=t2*t;t4=t3*t;si2=si*si;si3=si2*si;ri2=ri*ri;ri3=ri2*ri;
calc=A3*C*ri;
calc=calc + 3*A2*B*C*ri;
calc=calc + 3*A*B2*C*ri;
calc=calc + B3*C*ri;
calc=calc + 2*A2*C2*ri;
calc=calc + 4*A*B*C2*ri;
calc=calc + 2*B2*C2*ri;
calc=calc + A*C3*ri;
calc=calc + B*C3*ri;
calc=calc - A3*C*ri2;
calc=calc - 3*A2*B*C*ri2;
calc=calc - 3*A*B2*C*ri2;
calc=calc - B3*C*ri2;
calc=calc - 2*A2*C2*ri2;
calc=calc - 4*A*B*C2*ri2;
calc=calc - 2*B2*C2*ri2;
calc=calc - A2*B*C*si;
calc=calc - 2*A*B2*C*si;
calc=calc - B3*C*si;
calc=calc - 2*A*B*C2*si;
calc=calc - 2*B2*C2*si;
calc=calc - B*C3*si;
calc=calc - A3*C*ri*si;
calc=calc - A2*B*C*ri*si;
calc=calc + A*B2*C*ri*si;
calc=calc + B3*C*ri*si;
calc=calc - 2*A2*C2*ri*si;
calc=calc + 2*B2*C2*ri*si;
calc=calc - A*C3*ri*si;
calc=calc + A3*C*ri2*si;
calc=calc + 2*A2*B*C*ri2*si;
calc=calc + A*B2*C*ri2*si;
calc=calc + 2*A2*C2*ri2*si;
calc=calc + 2*A*B*C2*ri2*si;
calc=calc - A3*t;
calc=calc - 3*A2*B*t;
calc=calc - 3*A*B2*t;
calc=calc - B3*t;
calc=calc - 3*A2*C*t;
calc=calc - 6*A*B*C*t;
calc=calc - 3*B2*C*t;
calc=calc - 3*A*C2*t;
calc=calc - 3*B*C2*t;
calc=calc - C3*t;
calc=calc + A3*ri*t;
calc=calc + 3*A2*B*ri*t;
calc=calc + 3*A*B2*ri*t;
calc=calc + B3*ri*t;
calc=calc + 2*A2*C*ri*t;
calc=calc - A3*C*ri*t;
calc=calc + 4*A*B*C*ri*t;
calc=calc - 3*A2*B*C*ri*t;
calc=calc + 2*B2*C*ri*t;
calc=calc - 3*A*B2*C*ri*t;
calc=calc - B3*C*ri*t;
calc=calc + A*C2*ri*t;
calc=calc - 2*A2*C2*ri*t;
calc=calc + B*C2*ri*t;
calc=calc - 4*A*B*C2*ri*t;
calc=calc - 2*B2*C2*ri*t;
calc=calc - A*C3*ri*t;
calc=calc - B*C3*ri*t;
calc=calc + A2*C*ri2*t;
calc=calc + A3*C*ri2*t;
calc=calc + 2*A*B*C*ri2*t;
calc=calc + 3*A2*B*C*ri2*t;
calc=calc + B2*C*ri2*t;
calc=calc + 3*A*B2*C*ri2*t;
calc=calc + B3*C*ri2*t;
calc=calc + A*C2*ri2*t;
calc=calc + 3*A2*C2*ri2*t;
calc=calc + B*C2*ri2*t;
calc=calc + 4*A*B*C2*ri2*t;
calc=calc + B2*C2*ri2*t;
calc=calc + A3*si*t;
calc=calc + 2*A2*B*si*t;
calc=calc + A*B2*si*t;
calc=calc + 3*A2*C*si*t;
calc=calc + 4*A*B*C*si*t;
calc=calc + A2*B*C*si*t;
calc=calc + B2*C*si*t;
calc=calc + 2*A*B2*C*si*t;
calc=calc + B3*C*si*t;
calc=calc + 3*A*C2*si*t;
calc=calc + 2*B*C2*si*t;
calc=calc + 2*A*B*C2*si*t;
calc=calc + 2*B2*C2*si*t;
calc=calc + C3*si*t;
calc=calc + B*C3*si*t;
calc=calc - A3*ri*si*t;
calc=calc - A2*B*ri*si*t;
calc=calc + A*B2*ri*si*t;
calc=calc + B3*ri*si*t;
calc=calc - 2*A2*C*ri*si*t;
calc=calc + A3*C*ri*si*t;
calc=calc - A*B*C*ri*si*t;
calc=calc + 2*A2*B*C*ri*si*t;
calc=calc + B2*C*ri*si*t;
calc=calc - B3*C*ri*si*t;
calc=calc - A*C2*ri*si*t;
calc=calc + 2*A2*C2*ri*si*t;
calc=calc - B2*C2*ri*si*t;
calc=calc + A*C3*ri*si*t;
calc=calc - A2*B*ri2*si*t;
calc=calc - 2*A*B2*ri2*si*t;
calc=calc - B3*ri2*si*t;
calc=calc - A2*C*ri2*si*t;
calc=calc - A3*C*ri2*si*t;
calc=calc - 3*A*B*C*ri2*si*t;
calc=calc - 3*A2*B*C*ri2*si*t;
calc=calc - 2*B2*C*ri2*si*t;
calc=calc - 2*A*B2*C*ri2*si*t;
calc=calc - A*C2*ri2*si*t;
calc=calc - 3*A2*C2*ri2*si*t;
calc=calc - B*C2*ri2*si*t;
calc=calc - 2*A*B*C2*ri2*si*t;
calc=calc + A3*t2;
calc=calc + 3*A2*B*t2;
calc=calc + 3*A*B2*t2;
calc=calc + B3*t2;
calc=calc + 3*A2*C*t2;
calc=calc + 6*A*B*C*t2;
calc=calc + 3*B2*C*t2;
calc=calc + 3*A*C2*t2;
calc=calc + 3*B*C2*t2;
calc=calc + C3*t2;
calc=calc - A3*ri*t2;
calc=calc - 3*A2*B*ri*t2;
calc=calc - 3*A*B2*ri*t2;
calc=calc - B3*ri*t2;
calc=calc - 2*A2*C*ri*t2;
calc=calc - 4*A*B*C*ri*t2;
calc=calc - 2*B2*C*ri*t2;
calc=calc - A*C2*ri*t2;
calc=calc - B*C2*ri*t2;
calc=calc - A2*C*ri2*t2;
calc=calc - 2*A*B*C*ri2*t2;
calc=calc - B2*C*ri2*t2;
calc=calc - A*C2*ri2*t2;
calc=calc - A2*C2*ri2*t2;
calc=calc - B*C2*ri2*t2;
calc=calc - A3*si*t2;
calc=calc - 2*A2*B*si*t2;
calc=calc - A*B2*si*t2;
calc=calc - 3*A2*C*si*t2;
calc=calc - 4*A*B*C*si*t2;
calc=calc - B2*C*si*t2;
calc=calc - 3*A*C2*si*t2;
calc=calc - 2*B*C2*si*t2;
calc=calc - C3*si*t2;
calc=calc + A3*ri*si*t2;
calc=calc + A2*B*ri*si*t2;
calc=calc - A*B2*ri*si*t2;
calc=calc - B3*ri*si*t2;
calc=calc + 2*A2*C*ri*si*t2;
calc=calc + A*B*C*ri*si*t2;
calc=calc - A2*B*C*ri*si*t2;
calc=calc - B2*C*ri*si*t2;
calc=calc + A*C2*ri*si*t2;
calc=calc + A2*B*ri2*si*t2;
calc=calc + 2*A*B2*ri2*si*t2;
calc=calc + B3*ri2*si*t2;
calc=calc + A2*C*ri2*si*t2;
calc=calc + 3*A*B*C*ri2*si*t2;
calc=calc + A2*B*C*ri2*si*t2;
calc=calc + 2*B2*C*ri2*si*t2;
calc=calc + A*C2*ri2*si*t2;
calc=calc + A2*C2*ri2*si*t2;
calc=calc + B*C2*ri2*si*t2; return calc}

function crossterm(si, ri, A, B, C, t) {
A2=A*A;A3=A2*A;A4=A3*A;B2=B*B;B3=B2*B;B4=B3*B;C2=C*C;C3=C2*C;C4=C3*C;
t2=t*t;t3=t2*t;t4=t3*t;si2=si*si;si3=si2*si;ri2=ri*ri;ri3=ri2*ri;
calc= A2*B*C3*ri2*si;
calc=calc + 2*A*B2*C3*ri2*si;
calc=calc + B3*C3*ri2*si;
calc=calc - A2*B*C3*ri3*si;
calc=calc - 2*A*B2*C3*ri3*si;
calc=calc - B3*C3*ri3*si;
calc=calc - 2*A*B2*C3*ri*si2;
calc=calc - 2*B3*C3*ri*si2;
calc=calc - 2*A2*B*C3*ri2*si2;
calc=calc + 2*B3*C3*ri2*si2;
calc=calc + 2*A2*B*C3*ri3*si2;
calc=calc + 2*A*B2*C3*ri3*si2;
calc=calc + B3*C3*si3;
calc=calc + 2*A*B2*C3*ri*si3;
calc=calc - B3*C3*ri*si3;
calc=calc + A2*B*C3*ri2*si3;
calc=calc - 2*A*B2*C3*ri2*si3;
calc=calc - A2*B*C3*ri3*si3;
calc=calc - A4*C*ri*t;
calc=calc - 4*A3*B*C*ri*t;
calc=calc - 6*A2*B2*C*ri*t;
calc=calc - 4*A*B3*C*ri*t;
calc=calc - B4*C*ri*t;
calc=calc - 3*A3*C2*ri*t;
calc=calc - 9*A2*B*C2*ri*t;
calc=calc - 9*A*B2*C2*ri*t;
calc=calc - 3*B3*C2*ri*t;
calc=calc - 3*A2*C3*ri*t;
calc=calc - 6*A*B*C3*ri*t;
calc=calc - 3*B2*C3*ri*t;
calc=calc - A*C4*ri*t;
calc=calc - B*C4*ri*t;
calc=calc + 2*A4*C*ri2*t;
calc=calc + 8*A3*B*C*ri2*t;
calc=calc + 12*A2*B2*C*ri2*t;
calc=calc + 8*A*B3*C*ri2*t;
calc=calc + 2*B4*C*ri2*t;
calc=calc + 5*A3*C2*ri2*t;
calc=calc + 15*A2*B*C2*ri2*t;
calc=calc - A3*B*C2*ri2*t;
calc=calc + 15*A*B2*C2*ri2*t;
calc=calc - 3*A2*B2*C2*ri2*t;
calc=calc + 5*B3*C2*ri2*t;
calc=calc - 3*A*B3*C2*ri2*t;
calc=calc - B4*C2*ri2*t;
calc=calc + 3*A2*C3*ri2*t;
calc=calc + 6*A*B*C3*ri2*t;
calc=calc - 2*A2*B*C3*ri2*t;
calc=calc + 3*B2*C3*ri2*t;
calc=calc - 4*A*B2*C3*ri2*t;
calc=calc - 2*B3*C3*ri2*t;
calc=calc - A*B*C4*ri2*t;
calc=calc - B2*C4*ri2*t;
calc=calc - A4*C*ri3*t;
calc=calc - 4*A3*B*C*ri3*t;
calc=calc - 6*A2*B2*C*ri3*t;
calc=calc - 4*A*B3*C*ri3*t;
calc=calc - B4*C*ri3*t;
calc=calc - 2*A3*C2*ri3*t;
calc=calc - 6*A2*B*C2*ri3*t;
calc=calc + A3*B*C2*ri3*t;
calc=calc - 6*A*B2*C2*ri3*t;
calc=calc + 3*A2*B2*C2*ri3*t;
calc=calc - 2*B3*C2*ri3*t;
calc=calc + 3*A*B3*C2*ri3*t;
calc=calc + B4*C2*ri3*t;
calc=calc + 2*A2*B*C3*ri3*t;
calc=calc + 4*A*B2*C3*ri3*t;
calc=calc + 2*B3*C3*ri3*t;
calc=calc + A3*B*C*si*t;
calc=calc + 3*A2*B2*C*si*t;
calc=calc + 3*A*B3*C*si*t;
calc=calc + B4*C*si*t;
calc=calc + 3*A2*B*C2*si*t;
calc=calc + 6*A*B2*C2*si*t;
calc=calc + 3*B3*C2*si*t;
calc=calc + 3*A*B*C3*si*t;
calc=calc + 3*B2*C3*si*t;
calc=calc + B*C4*si*t;
calc=calc + 2*A4*C*ri*si*t;
calc=calc + 4*A3*B*C*ri*si*t;
calc=calc + A4*B*C*ri*si*t;
calc=calc + 3*A3*B2*C*ri*si*t;
calc=calc - 4*A*B3*C*ri*si*t;
calc=calc + 3*A2*B3*C*ri*si*t;
calc=calc - 2*B4*C*ri*si*t;
calc=calc + A*B4*C*ri*si*t;
calc=calc + 6*A3*C2*ri*si*t;
calc=calc + 7*A2*B*C2*ri*si*t;
calc=calc + 3*A3*B*C2*ri*si*t;
calc=calc - 4*A*B2*C2*ri*si*t;
calc=calc + 8*A2*B2*C2*ri*si*t;
calc=calc - 5*B3*C2*ri*si*t;
calc=calc + 7*A*B3*C2*ri*si*t;
calc=calc + 2*B4*C2*ri*si*t;
calc=calc + 6*A2*C3*ri*si*t;
calc=calc + 4*A*B*C3*ri*si*t;
calc=calc + 3*A2*B*C3*ri*si*t;
calc=calc - 2*B2*C3*ri*si*t;
calc=calc + 7*A*B2*C3*ri*si*t;
calc=calc + 4*B3*C3*ri*si*t;
calc=calc + 2*A*C4*ri*si*t;
calc=calc + B*C4*ri*si*t;
calc=calc + A*B*C4*ri*si*t;
calc=calc + 2*B2*C4*ri*si*t;
calc=calc - 4*A4*C*ri2*si*t;
calc=calc - 11*A3*B*C*ri2*si*t;
calc=calc - 2*A4*B*C*ri2*si*t;
calc=calc - 9*A2*B2*C*ri2*si*t;
calc=calc - 6*A3*B2*C*ri2*si*t;
calc=calc - A*B3*C*ri2*si*t;
calc=calc - 6*A2*B3*C*ri2*si*t;
calc=calc + B4*C*ri2*si*t;
calc=calc - 2*A*B4*C*ri2*si*t;
calc=calc - 9*A3*C2*ri2*si*t;
calc=calc - 16*A2*B*C2*ri2*si*t;
calc=calc - 2*A3*B*C2*ri2*si*t;
calc=calc - 5*A*B2*C2*ri2*si*t;
calc=calc - 6*A2*B2*C2*ri2*si*t;
calc=calc + 2*B3*C2*ri2*si*t;
calc=calc - 6*A*B3*C2*ri2*si*t;
calc=calc - 2*B4*C2*ri2*si*t;
calc=calc - 4*A2*C3*ri2*si*t;
calc=calc - 5*A*B*C3*ri2*si*t;
calc=calc - A2*B*C3*ri2*si*t;
calc=calc - B2*C3*ri2*si*t;
calc=calc - 5*A*B2*C3*ri2*si*t;
calc=calc - 4*B3*C3*ri2*si*t;
calc=calc + A*C4*ri2*si*t;
calc=calc + B*C4*ri2*si*t;
calc=calc + 2*A*B*C4*ri2*si*t;
calc=calc + B2*C4*ri2*si*t;
calc=calc + 2*A4*C*ri3*si*t;
calc=calc + 6*A3*B*C*ri3*si*t;
calc=calc + A4*B*C*ri3*si*t;
calc=calc + 6*A2*B2*C*ri3*si*t;
calc=calc + 3*A3*B2*C*ri3*si*t;
calc=calc + 2*A*B3*C*ri3*si*t;
calc=calc + 3*A2*B3*C*ri3*si*t;
calc=calc + A*B4*C*ri3*si*t;
calc=calc + 3*A3*C2*ri3*si*t;
calc=calc + 6*A2*B*C2*ri3*si*t;
calc=calc - A3*B*C2*ri3*si*t;
calc=calc + 3*A*B2*C2*ri3*si*t;
calc=calc - 2*A2*B2*C2*ri3*si*t;
calc=calc - A*B3*C2*ri3*si*t;
calc=calc - 2*A2*C3*ri3*si*t;
calc=calc - 2*A*B*C3*ri3*si*t;
calc=calc - 2*A2*B*C3*ri3*si*t;
calc=calc - 2*A*B2*C3*ri3*si*t;
calc=calc - A3*B*C*si2*t;
calc=calc - 2*A2*B2*C*si2*t;
calc=calc - A3*B2*C*si2*t;
calc=calc - A*B3*C*si2*t;
calc=calc - 2*A2*B3*C*si2*t;
calc=calc - A*B4*C*si2*t;
calc=calc - 3*A2*B*C2*si2*t;
calc=calc - 3*A*B2*C2*si2*t;
calc=calc - 3*A2*B2*C2*si2*t;
calc=calc - 4*A*B3*C2*si2*t;
calc=calc - B4*C2*si2*t;
calc=calc - 3*A*B*C3*si2*t;
calc=calc - B2*C3*si2*t;
calc=calc - 3*A*B2*C3*si2*t;
calc=calc - 2*B3*C3*si2*t;
calc=calc - B*C4*si2*t;
calc=calc - B2*C4*si2*t;
calc=calc - A4*C*ri*si2*t;
calc=calc + A3*B*C*ri*si2*t;
calc=calc - 2*A4*B*C*ri*si2*t;
calc=calc + 5*A2*B2*C*ri*si2*t;
calc=calc - 2*A3*B2*C*ri*si2*t;
calc=calc + 3*A*B3*C*ri*si2*t;
calc=calc + 2*A2*B3*C*ri*si2*t;
calc=calc + 2*A*B4*C*ri*si2*t;
calc=calc - 3*A3*C2*ri*si2*t;
calc=calc + 4*A2*B*C2*ri*si2*t;
calc=calc - 6*A3*B*C2*ri*si2*t;
calc=calc + 8*A*B2*C2*ri*si2*t;
calc=calc - 8*A2*B2*C2*ri*si2*t;
calc=calc + B3*C2*ri*si2*t;
calc=calc - A*B3*C2*ri*si2*t;
calc=calc + B4*C2*ri*si2*t;
calc=calc - 3*A2*C3*ri*si2*t;
calc=calc + 2*A*B*C3*ri*si2*t;
calc=calc - 6*A2*B*C3*ri*si2*t;
calc=calc + 3*B2*C3*ri*si2*t;
calc=calc - 2*A*B2*C3*ri*si2*t;
calc=calc + 2*B3*C3*ri*si2*t;
calc=calc - A*C4*ri*si2*t;
calc=calc - B*C4*ri*si2*t;
calc=calc - 2*A*B*C4*ri*si2*t;
calc=calc - 2*B2*C4*ri*si2*t;
calc=calc + 2*A4*C*ri2*si2*t;
calc=calc + A3*B*C*ri2*si2*t;
calc=calc + 4*A4*B*C*ri2*si2*t;
calc=calc - 4*A2*B2*C*ri2*si2*t;
calc=calc + 7*A3*B2*C*ri2*si2*t;
calc=calc - 3*A*B3*C*ri2*si2*t;
calc=calc + 2*A2*B3*C*ri2*si2*t;
calc=calc - A*B4*C*ri2*si2*t;
calc=calc + 3*A3*C2*ri2*si2*t;
calc=calc - 4*A2*B*C2*ri2*si2*t;
calc=calc + 7*A3*B*C2*ri2*si2*t;
calc=calc - 8*A*B2*C2*ri2*si2*t;
calc=calc + 14*A2*B2*C2*ri2*si2*t;
calc=calc - B3*C2*ri2*si2*t;
calc=calc + 7*A*B3*C2*ri2*si2*t;
calc=calc - A2*C3*ri2*si2*t;
calc=calc - A*B*C3*ri2*si2*t;
calc=calc + 8*A2*B*C3*ri2*si2*t;
calc=calc - 2*B2*C3*ri2*si2*t;
calc=calc + 7*A*B2*C3*ri2*si2*t;
calc=calc - 2*A*C4*ri2*si2*t;
calc=calc - B*C4*ri2*si2*t;
calc=calc - A*B*C4*ri2*si2*t;
calc=calc - A4*C*ri3*si2*t;
calc=calc - A3*B*C*ri3*si2*t;
calc=calc - 2*A4*B*C*ri3*si2*t;
calc=calc + A2*B2*C*ri3*si2*t;
calc=calc - 4*A3*B2*C*ri3*si2*t;
calc=calc + A*B3*C*ri3*si2*t;
calc=calc - 2*A2*B3*C*ri3*si2*t;
calc=calc + 3*A2*B*C2*ri3*si2*t;
calc=calc - A3*B*C2*ri3*si2*t;
calc=calc + 3*A*B2*C2*ri3*si2*t;
calc=calc - 3*A2*B2*C2*ri3*si2*t;
calc=calc - 2*A*B3*C2*ri3*si2*t;
calc=calc + 4*A2*C3*ri3*si2*t;
calc=calc + 2*A*B*C3*ri3*si2*t;
calc=calc - 2*A2*B*C3*ri3*si2*t;
calc=calc - 2*A*B2*C3*ri3*si2*t;
calc=calc - A2*B2*C*si3*t;
calc=calc + A3*B2*C*si3*t;
calc=calc - A*B3*C*si3*t;
calc=calc + A2*B3*C*si3*t;
calc=calc - 3*A*B2*C2*si3*t;
calc=calc + 3*A2*B2*C2*si3*t;
calc=calc - B3*C2*si3*t;
calc=calc + 4*A*B3*C2*si3*t;
calc=calc - 2*B2*C3*si3*t;
calc=calc + 3*A*B2*C3*si3*t;
calc=calc + B2*C4*si3*t;
calc=calc - A3*B*C*ri*si3*t;
calc=calc + A4*B*C*ri*si3*t;
calc=calc + A2*B2*C*ri*si3*t;
calc=calc - A3*B2*C*ri*si3*t;
calc=calc + 2*A*B3*C*ri*si3*t;
calc=calc - 2*A2*B3*C*ri*si3*t;
calc=calc - 2*A2*B*C2*ri*si3*t;
calc=calc + 3*A3*B*C2*ri*si3*t;
calc=calc + 5*A*B2*C2*ri*si3*t;
calc=calc + B3*C2*ri*si3*t;
calc=calc - 6*A*B3*C2*ri*si3*t;
calc=calc + 3*A2*B*C3*ri*si3*t;
calc=calc + 2*B2*C3*ri*si3*t;
calc=calc - 5*A*B2*C3*ri*si3*t;
calc=calc + B*C4*ri*si3*t;
calc=calc + A*B*C4*ri*si3*t;
calc=calc + 2*A3*B*C*ri2*si3*t;
calc=calc - 2*A4*B*C*ri2*si3*t;
calc=calc + A2*B2*C*ri2*si3*t;
calc=calc - A3*B2*C*ri2*si3*t;
calc=calc - A*B3*C*ri2*si3*t;
calc=calc + A2*B3*C*ri2*si3*t;
calc=calc + A3*C2*ri2*si3*t;
calc=calc + 5*A2*B*C2*ri2*si3*t;
calc=calc - 4*A3*B*C2*ri2*si3*t;
calc=calc - 2*A*B2*C2*ri2*si3*t;
calc=calc - 5*A2*B2*C2*ri2*si3*t;
calc=calc + 2*A*B3*C2*ri2*si3*t;
calc=calc + 2*A2*C3*ri2*si3*t;
calc=calc - 5*A2*B*C3*ri2*si3*t;
calc=calc + 2*A*B2*C3*ri2*si3*t;
calc=calc + A*C4*ri2*si3*t;
calc=calc - A3*B*C*ri3*si3*t;
calc=calc + A4*B*C*ri3*si3*t;
calc=calc - A2*B2*C*ri3*si3*t;
calc=calc + A3*B2*C*ri3*si3*t;
calc=calc - A3*C2*ri3*si3*t;
calc=calc - 3*A2*B*C2*ri3*si3*t;
calc=calc + A3*B*C2*ri3*si3*t;
calc=calc + 2*A2*B2*C2*ri3*si3*t;
calc=calc - 2*A2*C3*ri3*si3*t;
calc=calc + 2*A2*B*C3*ri3*si3*t;
calc=calc + A4*t2;
calc=calc + 4*A3*B*t2;
calc=calc + 6*A2*B2*t2;
calc=calc + 4*A*B3*t2;
calc=calc + B4*t2;
calc=calc + 4*A3*C*t2;
calc=calc + 12*A2*B*C*t2;
calc=calc + 12*A*B2*C*t2;
calc=calc + 4*B3*C*t2;
calc=calc + 6*A2*C2*t2;
calc=calc + 12*A*B*C2*t2;
calc=calc + 6*B2*C2*t2;
calc=calc + 4*A*C3*t2;
calc=calc + 4*B*C3*t2;
calc=calc + C4*t2;
calc=calc - 2*A4*ri*t2;
calc=calc - 8*A3*B*ri*t2;
calc=calc - 12*A2*B2*ri*t2;
calc=calc - 8*A*B3*ri*t2;
calc=calc - 2*B4*ri*t2;
calc=calc - 6*A3*C*ri*t2;
calc=calc + 2*A4*C*ri*t2;
calc=calc - 18*A2*B*C*ri*t2;
calc=calc + 9*A3*B*C*ri*t2;
calc=calc - 18*A*B2*C*ri*t2;
calc=calc + 15*A2*B2*C*ri*t2;
calc=calc - 6*B3*C*ri*t2;
calc=calc + 11*A*B3*C*ri*t2;
calc=calc + 3*B4*C*ri*t2;
calc=calc - 6*A2*C2*ri*t2;
calc=calc + 6*A3*C2*ri*t2;
calc=calc - 12*A*B*C2*ri*t2;
calc=calc + 21*A2*B*C2*ri*t2;
calc=calc - 6*B2*C2*ri*t2;
calc=calc + 24*A*B2*C2*ri*t2;
calc=calc + 9*B3*C2*ri*t2;
calc=calc - 2*A*C3*ri*t2;
calc=calc + 6*A2*C3*ri*t2;
calc=calc - 2*B*C3*ri*t2;
calc=calc + 15*A*B*C3*ri*t2;
calc=calc + 9*B2*C3*ri*t2;
calc=calc + 2*A*C4*ri*t2;
calc=calc + 3*B*C4*ri*t2;
calc=calc + A4*ri2*t2;
calc=calc + 4*A3*B*ri2*t2;
calc=calc + 6*A2*B2*ri2*t2;
calc=calc + 4*A*B3*ri2*t2;
calc=calc + B4*ri2*t2;
calc=calc + A3*C*ri2*t2;
calc=calc - 4*A4*C*ri2*t2;
calc=calc + 3*A2*B*C*ri2*t2;
calc=calc - 17*A3*B*C*ri2*t2;
calc=calc + 3*A*B2*C*ri2*t2;
calc=calc - 27*A2*B2*C*ri2*t2;
calc=calc + B3*C*ri2*t2;
calc=calc - 19*A*B3*C*ri2*t2;
calc=calc - 5*B4*C*ri2*t2;
calc=calc - A2*C2*ri2*t2;
calc=calc - 11*A3*C2*ri2*t2;
calc=calc - 2*A*B*C2*ri2*t2;
calc=calc - 33*A2*B*C2*ri2*t2;
calc=calc + A3*B*C2*ri2*t2;
calc=calc - B2*C2*ri2*t2;
calc=calc - 33*A*B2*C2*ri2*t2;
calc=calc + 3*A2*B2*C2*ri2*t2;
calc=calc - 11*B3*C2*ri2*t2;
calc=calc + 3*A*B3*C2*ri2*t2;
calc=calc + B4*C2*ri2*t2;
calc=calc - A*C3*ri2*t2;
calc=calc - 7*A2*C3*ri2*t2;
calc=calc - B*C3*ri2*t2;
calc=calc - 13*A*B*C3*ri2*t2;
calc=calc + 2*A2*B*C3*ri2*t2;
calc=calc - 6*B2*C3*ri2*t2;
calc=calc + 4*A*B2*C3*ri2*t2;
calc=calc + 2*B3*C3*ri2*t2;
calc=calc + A*B*C4*ri2*t2;
calc=calc + B2*C4*ri2*t2;
calc=calc + A3*C*ri3*t2;
calc=calc + 2*A4*C*ri3*t2;
calc=calc + 3*A2*B*C*ri3*t2;
calc=calc + 8*A3*B*C*ri3*t2;
calc=calc + 3*A*B2*C*ri3*t2;
calc=calc + 12*A2*B2*C*ri3*t2;
calc=calc + B3*C*ri3*t2;
calc=calc + 8*A*B3*C*ri3*t2;
calc=calc + 2*B4*C*ri3*t2;
calc=calc + A2*C2*ri3*t2;
calc=calc + 5*A3*C2*ri3*t2;
calc=calc + 2*A*B*C2*ri3*t2;
calc=calc + 12*A2*B*C2*ri3*t2;
calc=calc - A3*B*C2*ri3*t2;
calc=calc + B2*C2*ri3*t2;
calc=calc + 9*A*B2*C2*ri3*t2;
calc=calc - 3*A2*B2*C2*ri3*t2;
calc=calc + 2*B3*C2*ri3*t2;
calc=calc - 3*A*B3*C2*ri3*t2;
calc=calc - B4*C2*ri3*t2;
calc=calc - A*B*C3*ri3*t2;
calc=calc - 3*A2*B*C3*ri3*t2;
calc=calc - B2*C3*ri3*t2;
calc=calc - 4*A*B2*C3*ri3*t2;
calc=calc - B3*C3*ri3*t2;
calc=calc - 2*A4*si*t2;
calc=calc - 6*A3*B*si*t2;
calc=calc - A4*B*si*t2;
calc=calc - 6*A2*B2*si*t2;
calc=calc - 3*A3*B2*si*t2;
calc=calc - 2*A*B3*si*t2;
calc=calc - 3*A2*B3*si*t2;
calc=calc - A*B4*si*t2;
calc=calc - 8*A3*C*si*t2;
calc=calc - 18*A2*B*C*si*t2;
calc=calc - 6*A3*B*C*si*t2;
calc=calc - 12*A*B2*C*si*t2;
calc=calc - 15*A2*B2*C*si*t2;
calc=calc - 2*B3*C*si*t2;
calc=calc - 12*A*B3*C*si*t2;
calc=calc - 3*B4*C*si*t2;
calc=calc - 12*A2*C2*si*t2;
calc=calc - 18*A*B*C2*si*t2;
calc=calc - 12*A2*B*C2*si*t2;
calc=calc - 6*B2*C2*si*t2;
calc=calc - 21*A*B2*C2*si*t2;
calc=calc - 9*B3*C2*si*t2;
calc=calc - 8*A*C3*si*t2;
calc=calc - 6*B*C3*si*t2;
calc=calc - 10*A*B*C3*si*t2;
calc=calc - 9*B2*C3*si*t2;
calc=calc - 2*C4*si*t2;
calc=calc - 3*B*C4*si*t2;
calc=calc + 4*A4*ri*si*t2;
calc=calc + 11*A3*B*ri*si*t2;
calc=calc + 2*A4*B*ri*si*t2;
calc=calc + 9*A2*B2*ri*si*t2;
calc=calc + 6*A3*B2*ri*si*t2;
calc=calc + A*B3*ri*si*t2;
calc=calc + 6*A2*B3*ri*si*t2;
calc=calc - B4*ri*si*t2;
calc=calc + 2*A*B4*ri*si*t2;
calc=calc + 11*A3*C*ri*si*t2;
calc=calc - 4*A4*C*ri*si*t2;
calc=calc + 22*A2*B*C*ri*si*t2;
calc=calc - 6*A3*B*C*ri*si*t2;
calc=calc - A4*B*C*ri*si*t2;
calc=calc + 11*A*B2*C*ri*si*t2;
calc=calc + 5*A2*B2*C*ri*si*t2;
calc=calc - 3*A3*B2*C*ri*si*t2;
calc=calc + 12*A*B3*C*ri*si*t2;
calc=calc - 3*A2*B3*C*ri*si*t2;
calc=calc + 5*B4*C*ri*si*t2;
calc=calc - A*B4*C*ri*si*t2;
calc=calc + 9*A2*C2*ri*si*t2;
calc=calc - 12*A3*C2*ri*si*t2;
calc=calc + 11*A*B*C2*ri*si*t2;
calc=calc - 16*A2*B*C2*ri*si*t2;
calc=calc - 3*A3*B*C2*ri*si*t2;
calc=calc + 2*B2*C2*ri*si*t2;
calc=calc + 4*A*B2*C2*ri*si*t2;
calc=calc - 8*A2*B2*C2*ri*si*t2;
calc=calc + 8*B3*C2*ri*si*t2;
calc=calc - 7*A*B3*C2*ri*si*t2;
calc=calc - 2*B4*C2*ri*si*t2;
calc=calc + A*C3*ri*si*t2;
calc=calc - 12*A2*C3*ri*si*t2;
calc=calc - 12*A*B*C3*ri*si*t2;
calc=calc - 3*A2*B*C3*ri*si*t2;
calc=calc - B2*C3*ri*si*t2;
calc=calc - 7*A*B2*C3*ri*si*t2;
calc=calc - 4*B3*C3*ri*si*t2;
calc=calc - C4*ri*si*t2;
calc=calc - 4*A*C4*ri*si*t2;
calc=calc - 4*B*C4*ri*si*t2;
calc=calc - A*B*C4*ri*si*t2;
calc=calc - 2*B2*C4*ri*si*t2;
calc=calc - 2*A4*ri2*si*t2;
calc=calc - 4*A3*B*ri2*si*t2;
calc=calc - A4*B*ri2*si*t2;
calc=calc - 3*A3*B2*ri2*si*t2;
calc=calc + 4*A*B3*ri2*si*t2;
calc=calc - 3*A2*B3*ri2*si*t2;
calc=calc + 2*B4*ri2*si*t2;
calc=calc - A*B4*ri2*si*t2;
calc=calc - A3*C*ri2*si*t2;
calc=calc + 8*A4*C*ri2*si*t2;
calc=calc + 2*A2*B*C*ri2*si*t2;
calc=calc + 26*A3*B*C*ri2*si*t2;
calc=calc + 2*A4*B*C*ri2*si*t2;
calc=calc + 7*A*B2*C*ri2*si*t2;
calc=calc + 25*A2*B2*C*ri2*si*t2;
calc=calc + 6*A3*B2*C*ri2*si*t2;
calc=calc + 4*B3*C*ri2*si*t2;
calc=calc + 4*A*B3*C*ri2*si*t2;
calc=calc + 6*A2*B3*C*ri2*si*t2;
calc=calc - 3*B4*C*ri2*si*t2;
calc=calc + 2*A*B4*C*ri2*si*t2;
calc=calc + 4*A2*C2*ri2*si*t2;
calc=calc + 20*A3*C2*ri2*si*t2;
calc=calc + 9*A*B*C2*ri2*si*t2;
calc=calc + 39*A2*B*C2*ri2*si*t2;
calc=calc + 2*A3*B*C2*ri2*si*t2;
calc=calc + 5*B2*C2*ri2*si*t2;
calc=calc + 18*A*B2*C2*ri2*si*t2;
calc=calc + 4*A2*B2*C2*ri2*si*t2;
calc=calc - B3*C2*ri2*si*t2;
calc=calc + 4*A*B3*C2*ri2*si*t2;
calc=calc + 2*B4*C2*ri2*si*t2;
calc=calc + 3*A*C3*ri2*si*t2;
calc=calc + 10*A2*C3*ri2*si*t2;
calc=calc + 3*B*C3*ri2*si*t2;
calc=calc + 13*A*B*C3*ri2*si*t2;
calc=calc + A2*B*C3*ri2*si*t2;
calc=calc + 3*B2*C3*ri2*si*t2;
calc=calc + 3*A*B2*C3*ri2*si*t2;
calc=calc + B3*C3*ri2*si*t2;
calc=calc - 2*A*C4*ri2*si*t2;
calc=calc - 2*B*C4*ri2*si*t2;
calc=calc - 2*A*B*C4*ri2*si*t2;
calc=calc - B2*C4*ri2*si*t2;
calc=calc - A3*B*ri3*si*t2;
calc=calc - 3*A2*B2*ri3*si*t2;
calc=calc - 3*A*B3*ri3*si*t2;
calc=calc - B4*ri3*si*t2;
calc=calc - 2*A3*C*ri3*si*t2;
calc=calc - 4*A4*C*ri3*si*t2;
calc=calc - 6*A2*B*C*ri3*si*t2;
calc=calc - 14*A3*B*C*ri3*si*t2;
calc=calc - A4*B*C*ri3*si*t2;
calc=calc - 6*A*B2*C*ri3*si*t2;
calc=calc - 15*A2*B2*C*ri3*si*t2;
calc=calc - 3*A3*B2*C*ri3*si*t2;
calc=calc - 2*B3*C*ri3*si*t2;
calc=calc - 4*A*B3*C*ri3*si*t2;
calc=calc - 3*A2*B3*C*ri3*si*t2;
calc=calc + B4*C*ri3*si*t2;
calc=calc - A*B4*C*ri3*si*t2;
calc=calc - A2*C2*ri3*si*t2;
calc=calc - 8*A3*C2*ri3*si*t2;
calc=calc - 2*A*B*C2*ri3*si*t2;
calc=calc - 11*A2*B*C2*ri3*si*t2;
calc=calc + A3*B*C2*ri3*si*t2;
calc=calc - B2*C2*ri3*si*t2;
calc=calc - A*B2*C2*ri3*si*t2;
calc=calc + 4*A2*B2*C2*ri3*si*t2;
calc=calc + 2*B3*C2*ri3*si*t2;
calc=calc + 3*A*B3*C2*ri3*si*t2;
calc=calc + A*C3*ri3*si*t2;
calc=calc + 5*A2*C3*ri3*si*t2;
calc=calc + 6*A*B*C3*ri3*si*t2;
calc=calc + 5*A2*B*C3*ri3*si*t2;
calc=calc + B2*C3*ri3*si*t2;
calc=calc + 4*A*B2*C3*ri3*si*t2;
calc=calc + A4*si2*t2;
calc=calc + A3*B*si2*t2;
calc=calc + 2*A4*B*si2*t2;
calc=calc - A2*B2*si2*t2;
calc=calc + 4*A3*B2*si2*t2;
calc=calc - A*B3*si2*t2;
calc=calc + 2*A2*B3*si2*t2;
calc=calc + 4*A3*C*si2*t2;
calc=calc + 3*A2*B*C*si2*t2;
calc=calc + 10*A3*B*C*si2*t2;
calc=calc - 2*A*B2*C*si2*t2;
calc=calc + 18*A2*B2*C*si2*t2;
calc=calc + A3*B2*C*si2*t2;
calc=calc - B3*C*si2*t2;
calc=calc + 8*A*B3*C*si2*t2;
calc=calc + 2*A2*B3*C*si2*t2;
calc=calc + A*B4*C*si2*t2;
calc=calc + 6*A2*C2*si2*t2;
calc=calc + 3*A*B*C2*si2*t2;
calc=calc + 18*A2*B*C2*si2*t2;
calc=calc - B2*C2*si2*t2;
calc=calc + 21*A*B2*C2*si2*t2;
calc=calc + 3*A2*B2*C2*si2*t2;
calc=calc + 3*B3*C2*si2*t2;
calc=calc + 4*A*B3*C2*si2*t2;
calc=calc + B4*C2*si2*t2;
calc=calc + 4*A*C3*si2*t2;
calc=calc + B*C3*si2*t2;
calc=calc + 14*A*B*C3*si2*t2;
calc=calc + 7*B2*C3*si2*t2;
calc=calc + 3*A*B2*C3*si2*t2;
calc=calc + 2*B3*C3*si2*t2;
calc=calc + C4*si2*t2;
calc=calc + 4*B*C4*si2*t2;
calc=calc + B2*C4*si2*t2;
calc=calc - 2*A4*ri*si2*t2;
calc=calc - A3*B*ri*si2*t2;
calc=calc - 4*A4*B*ri*si2*t2;
calc=calc + 4*A2*B2*ri*si2*t2;
calc=calc - 7*A3*B2*ri*si2*t2;
calc=calc + 3*A*B3*ri*si2*t2;
calc=calc - 2*A2*B3*ri*si2*t2;
calc=calc + A*B4*ri*si2*t2;
calc=calc - 4*A3*C*ri*si2*t2;
calc=calc + 2*A4*C*ri*si2*t2;
calc=calc + 2*A2*B*C*ri*si2*t2;
calc=calc - 10*A3*B*C*ri*si2*t2;
calc=calc + 2*A4*B*C*ri*si2*t2;
calc=calc + 9*A*B2*C*ri*si2*t2;
calc=calc - 24*A2*B2*C*ri*si2*t2;
calc=calc + 2*A3*B2*C*ri*si2*t2;
calc=calc + 3*B3*C*ri*si2*t2;
calc=calc - 11*A*B3*C*ri*si2*t2;
calc=calc - 2*A2*B3*C*ri*si2*t2;
calc=calc + B4*C*ri*si2*t2;
calc=calc - 2*A*B4*C*ri*si2*t2;
calc=calc + 6*A3*C2*ri*si2*t2;
calc=calc + 7*A*B*C2*ri*si2*t2;
calc=calc - 14*A2*B*C2*ri*si2*t2;
calc=calc + 6*A3*B*C2*ri*si2*t2;
calc=calc + 5*B2*C2*ri*si2*t2;
calc=calc - 21*A*B2*C2*ri*si2*t2;
calc=calc + 10*A2*B2*C2*ri*si2*t2;
calc=calc - 3*B3*C2*ri*si2*t2;
calc=calc + 3*A*B3*C2*ri*si2*t2;
calc=calc - B4*C2*ri*si2*t2;
calc=calc + 4*A*C3*ri*si2*t2;
calc=calc + 6*A2*C3*ri*si2*t2;
calc=calc + 4*B*C3*ri*si2*t2;
calc=calc - 5*A*B*C3*ri*si2*t2;
calc=calc + 6*A2*B*C3*ri*si2*t2;
calc=calc - 4*B2*C3*ri*si2*t2;
calc=calc + 4*A*B2*C3*ri*si2*t2;
calc=calc + B3*C3*ri*si2*t2;
calc=calc + 2*C4*ri*si2*t2;
calc=calc + 2*A*C4*ri*si2*t2;
calc=calc + 3*B*C4*ri*si2*t2;
calc=calc + 2*A*B*C4*ri*si2*t2;
calc=calc + 2*B2*C4*ri*si2*t2;
calc=calc + A4*ri2*si2*t2;
calc=calc - A3*B*ri2*si2*t2;
calc=calc + 2*A4*B*ri2*si2*t2;
calc=calc - 5*A2*B2*ri2*si2*t2;
calc=calc + 2*A3*B2*ri2*si2*t2;
calc=calc - 3*A*B3*ri2*si2*t2;
calc=calc - 2*A2*B3*ri2*si2*t2;
calc=calc - 2*A*B4*ri2*si2*t2;
calc=calc - A3*C*ri2*si2*t2;
calc=calc - 4*A4*C*ri2*si2*t2;
calc=calc - 6*A2*B*C*ri2*si2*t2;
calc=calc - 5*A3*B*C*ri2*si2*t2;
calc=calc - 4*A4*B*C*ri2*si2*t2;
calc=calc - 7*A*B2*C*ri2*si2*t2;
calc=calc + 5*A2*B2*C*ri2*si2*t2;
calc=calc - 7*A3*B2*C*ri2*si2*t2;
calc=calc - 2*B3*C*ri2*si2*t2;
calc=calc + 5*A*B3*C*ri2*si2*t2;
calc=calc - 2*A2*B3*C*ri2*si2*t2;
calc=calc - B4*C*ri2*si2*t2;
calc=calc + A*B4*C*ri2*si2*t2;
calc=calc - 5*A2*C2*ri2*si2*t2;
calc=calc - 7*A3*C2*ri2*si2*t2;
calc=calc - 7*A*B*C2*ri2*si2*t2;
calc=calc + 5*A2*B*C2*ri2*si2*t2;
calc=calc - 7*A3*B*C2*ri2*si2*t2;
calc=calc - 4*B2*C2*ri2*si2*t2;
calc=calc + 8*A*B2*C2*ri2*si2*t2;
calc=calc - 14*A2*B2*C2*ri2*si2*t2;
calc=calc - 7*A*B3*C2*ri2*si2*t2;
calc=calc - 3*A*C3*ri2*si2*t2;
calc=calc + A2*C3*ri2*si2*t2;
calc=calc - 2*B*C3*ri2*si2*t2;
calc=calc - A*B*C3*ri2*si2*t2;
calc=calc - 8*A2*B*C3*ri2*si2*t2;
calc=calc + 3*B2*C3*ri2*si2*t2;
calc=calc - 7*A*B2*C3*ri2*si2*t2;
calc=calc + 4*A*C4*ri2*si2*t2;
calc=calc + 2*B*C4*ri2*si2*t2;
calc=calc + A*B*C4*ri2*si2*t2;
calc=calc + A3*B*ri3*si2*t2;
calc=calc + 2*A2*B2*ri3*si2*t2;
calc=calc + A3*B2*ri3*si2*t2;
calc=calc + A*B3*ri3*si2*t2;
calc=calc + 2*A2*B3*ri3*si2*t2;
calc=calc + A*B4*ri3*si2*t2;
calc=calc + A3*C*ri3*si2*t2;
calc=calc + 2*A4*C*ri3*si2*t2;
calc=calc + A2*B*C*ri3*si2*t2;
calc=calc + 5*A3*B*C*ri3*si2*t2;
calc=calc + 2*A4*B*C*ri3*si2*t2;
calc=calc + A2*B2*C*ri3*si2*t2;
calc=calc + 4*A3*B2*C*ri3*si2*t2;
calc=calc - 2*A*B3*C*ri3*si2*t2;
calc=calc + 2*A2*B3*C*ri3*si2*t2;
calc=calc - A2*C2*ri3*si2*t2;
calc=calc + A3*C2*ri3*si2*t2;
calc=calc - 3*A*B*C2*ri3*si2*t2;
calc=calc - 9*A2*B*C2*ri3*si2*t2;
calc=calc + A3*B*C2*ri3*si2*t2;
calc=calc - 8*A*B2*C2*ri3*si2*t2;
calc=calc + A2*B2*C2*ri3*si2*t2;
calc=calc - 2*A*C3*ri3*si2*t2;
calc=calc - 10*A2*C3*ri3*si2*t2;
calc=calc - 5*A*B*C3*ri3*si2*t2;
calc=calc - A2*B*C3*ri3*si2*t2;
calc=calc + A3*B*si3*t2;
calc=calc - A4*B*si3*t2;
calc=calc + A2*B2*si3*t2;
calc=calc - A3*B2*si3*t2;
calc=calc + 3*A2*B*C*si3*t2;
calc=calc - 4*A3*B*C*si3*t2;
calc=calc + 2*A*B2*C*si3*t2;
calc=calc - 3*A2*B2*C*si3*t2;
calc=calc - A3*B2*C*si3*t2;
calc=calc + 2*A*B3*C*si3*t2;
calc=calc + 3*A*B*C2*si3*t2;
calc=calc - 6*A2*B*C2*si3*t2;
calc=calc + B2*C2*si3*t2;
calc=calc - 3*A2*B2*C2*si3*t2;
calc=calc + 2*B3*C2*si3*t2;
calc=calc - 4*A*B3*C2*si3*t2;
calc=calc + B*C3*si3*t2;
calc=calc - 4*A*B*C3*si3*t2;
calc=calc + 2*B2*C3*si3*t2;
calc=calc - 3*A*B2*C3*si3*t2;
calc=calc - B3*C3*si3*t2;
calc=calc - B*C4*si3*t2;
calc=calc - B2*C4*si3*t2;
calc=calc - 2*A3*B*ri*si3*t2;
calc=calc + 2*A4*B*ri*si3*t2;
calc=calc - A2*B2*ri*si3*t2;
calc=calc + A3*B2*ri*si3*t2;
calc=calc + A*B3*ri*si3*t2;
calc=calc - A2*B3*ri*si3*t2;
calc=calc - A3*C*ri*si3*t2;
calc=calc - 6*A2*B*C*ri*si3*t2;
calc=calc + 7*A3*B*C*ri*si3*t2;
calc=calc - A4*B*C*ri*si3*t2;
calc=calc - 2*A*B2*C*ri*si3*t2;
calc=calc + 4*A2*B2*C*ri*si3*t2;
calc=calc + A3*B2*C*ri*si3*t2;
calc=calc - 6*A*B3*C*ri*si3*t2;
calc=calc - A2*B3*C*ri*si3*t2;
calc=calc - 3*A2*C2*ri*si3*t2;
calc=calc - 6*A*B*C2*ri*si3*t2;
calc=calc + 9*A2*B*C2*ri*si3*t2;
calc=calc - 3*A3*B*C2*ri*si3*t2;
calc=calc - B2*C2*ri*si3*t2;
calc=calc - 7*A*B2*C2*ri*si3*t2;
calc=calc - 2*A2*B2*C2*ri*si3*t2;
calc=calc - 2*B3*C2*ri*si3*t2;
calc=calc + 4*A*B3*C2*ri*si3*t2;
calc=calc - 3*A*C3*ri*si3*t2;
calc=calc - 2*B*C3*ri*si3*t2;
calc=calc + 2*A*B*C3*ri*si3*t2;
calc=calc - 3*A2*B*C3*ri*si3*t2;
calc=calc - 4*B2*C3*ri*si3*t2;
calc=calc + 3*A*B2*C3*ri*si3*t2;
calc=calc - C4*ri*si3*t2;
calc=calc - 2*B*C4*ri*si3*t2;
calc=calc - A*B*C4*ri*si3*t2;
calc=calc + A3*B*ri2*si3*t2;
calc=calc - A4*B*ri2*si3*t2;
calc=calc - A2*B2*ri2*si3*t2;
calc=calc + A3*B2*ri2*si3*t2;
calc=calc - 2*A*B3*ri2*si3*t2;
calc=calc + 2*A2*B3*ri2*si3*t2;
calc=calc + A3*C*ri2*si3*t2;
calc=calc + A2*B*C*ri2*si3*t2;
calc=calc - 4*A3*B*C*ri2*si3*t2;
calc=calc + 2*A4*B*C*ri2*si3*t2;
calc=calc - 3*A*B2*C*ri2*si3*t2;
calc=calc - 3*A2*B2*C*ri2*si3*t2;
calc=calc + A3*B2*C*ri2*si3*t2;
calc=calc + 4*A*B3*C*ri2*si3*t2;
calc=calc + 2*A2*B3*C*ri2*si3*t2;
calc=calc + 2*A2*C2*ri2*si3*t2;
calc=calc - 2*A3*C2*ri2*si3*t2;
calc=calc - 11*A2*B*C2*ri2*si3*t2;
calc=calc + 4*A3*B*C2*ri2*si3*t2;
calc=calc + 7*A*B2*C2*ri2*si3*t2;
calc=calc + 7*A2*B2*C2*ri2*si3*t2;
calc=calc + A*C3*ri2*si3*t2;
calc=calc - 4*A2*C3*ri2*si3*t2;
calc=calc + A*B*C3*ri2*si3*t2;
calc=calc + 5*A2*B*C3*ri2*si3*t2;
calc=calc - 2*A*C4*ri2*si3*t2;
calc=calc + A2*B2*ri3*si3*t2;
calc=calc - A3*B2*ri3*si3*t2;
calc=calc + A*B3*ri3*si3*t2;
calc=calc - A2*B3*ri3*si3*t2;
calc=calc + 2*A2*B*C*ri3*si3*t2;
calc=calc + A3*B*C*ri3*si3*t2;
calc=calc - A4*B*C*ri3*si3*t2;
calc=calc + 3*A*B2*C*ri3*si3*t2;
calc=calc + 2*A2*B2*C*ri3*si3*t2;
calc=calc - A3*B2*C*ri3*si3*t2;
calc=calc - A2*B3*C*ri3*si3*t2;
calc=calc + A2*C2*ri3*si3*t2;
calc=calc + 2*A3*C2*ri3*si3*t2;
calc=calc + 3*A*B*C2*ri3*si3*t2;
calc=calc + 8*A2*B*C2*ri3*si3*t2;
calc=calc - A3*B*C2*ri3*si3*t2;
calc=calc - 2*A2*B2*C2*ri3*si3*t2;
calc=calc + A*C3*ri3*si3*t2;
calc=calc + 5*A2*C3*ri3*si3*t2;
calc=calc - A2*B*C3*ri3*si3*t2;
calc=calc - 2*A4*t3;
calc=calc - 8*A3*B*t3;
calc=calc - 12*A2*B2*t3;
calc=calc - 8*A*B3*t3;
calc=calc - 2*B4*t3;
calc=calc - 8*A3*C*t3;
calc=calc - 24*A2*B*C*t3;
calc=calc - 24*A*B2*C*t3;
calc=calc - 8*B3*C*t3;
calc=calc - 12*A2*C2*t3;
calc=calc - 24*A*B*C2*t3;
calc=calc - 12*B2*C2*t3;
calc=calc - 8*A*C3*t3;
calc=calc - 8*B*C3*t3;
calc=calc - 2*C4*t3;
calc=calc + 4*A4*ri*t3;
calc=calc + 16*A3*B*ri*t3;
calc=calc + 24*A2*B2*ri*t3;
calc=calc + 16*A*B3*ri*t3;
calc=calc + 4*B4*ri*t3;
calc=calc + 12*A3*C*ri*t3;
calc=calc - A4*C*ri*t3;
calc=calc + 36*A2*B*C*ri*t3;
calc=calc - 5*A3*B*C*ri*t3;
calc=calc + 36*A*B2*C*ri*t3;
calc=calc - 9*A2*B2*C*ri*t3;
calc=calc + 12*B3*C*ri*t3;
calc=calc - 7*A*B3*C*ri*t3;
calc=calc - 2*B4*C*ri*t3;
calc=calc + 12*A2*C2*ri*t3;
calc=calc - 3*A3*C2*ri*t3;
calc=calc + 24*A*B*C2*ri*t3;
calc=calc - 12*A2*B*C2*ri*t3;
calc=calc + 12*B2*C2*ri*t3;
calc=calc - 15*A*B2*C2*ri*t3;
calc=calc - 6*B3*C2*ri*t3;
calc=calc + 4*A*C3*ri*t3;
calc=calc - 3*A2*C3*ri*t3;
calc=calc + 4*B*C3*ri*t3;
calc=calc - 9*A*B*C3*ri*t3;
calc=calc - 6*B2*C3*ri*t3;
calc=calc - A*C4*ri*t3;
calc=calc - 2*B*C4*ri*t3;
calc=calc - 2*A4*ri2*t3;
calc=calc - 8*A3*B*ri2*t3;
calc=calc - 12*A2*B2*ri2*t3;
calc=calc - 8*A*B3*ri2*t3;
calc=calc - 2*B4*ri2*t3;
calc=calc - 2*A3*C*ri2*t3;
calc=calc + 2*A4*C*ri2*t3;
calc=calc - 6*A2*B*C*ri2*t3;
calc=calc + 9*A3*B*C*ri2*t3;
calc=calc - 6*A*B2*C*ri2*t3;
calc=calc + 15*A2*B2*C*ri2*t3;
calc=calc - 2*B3*C*ri2*t3;
calc=calc + 11*A*B3*C*ri2*t3;
calc=calc + 3*B4*C*ri2*t3;
calc=calc + 2*A2*C2*ri2*t3;
calc=calc + 7*A3*C2*ri2*t3;
calc=calc + 4*A*B*C2*ri2*t3;
calc=calc + 19*A2*B*C2*ri2*t3;
calc=calc + 2*B2*C2*ri2*t3;
calc=calc + 18*A*B2*C2*ri2*t3;
calc=calc + 6*B3*C2*ri2*t3;
calc=calc + 2*A*C3*ri2*t3;
calc=calc + 5*A2*C3*ri2*t3;
calc=calc + 2*B*C3*ri2*t3;
calc=calc + 7*A*B*C3*ri2*t3;
calc=calc + 3*B2*C3*ri2*t3;
calc=calc - 2*A3*C*ri3*t3;
calc=calc - A4*C*ri3*t3;
calc=calc - 6*A2*B*C*ri3*t3;
calc=calc - 4*A3*B*C*ri3*t3;
calc=calc - 6*A*B2*C*ri3*t3;
calc=calc - 6*A2*B2*C*ri3*t3;
calc=calc - 2*B3*C*ri3*t3;
calc=calc - 4*A*B3*C*ri3*t3;
calc=calc - B4*C*ri3*t3;
calc=calc - 2*A2*C2*ri3*t3;
calc=calc - 4*A3*C2*ri3*t3;
calc=calc - 4*A*B*C2*ri3*t3;
calc=calc - 7*A2*B*C2*ri3*t3;
calc=calc - 2*B2*C2*ri3*t3;
calc=calc - 3*A*B2*C2*ri3*t3;
calc=calc + A*B*C3*ri3*t3;
calc=calc + A2*B*C3*ri3*t3;
calc=calc + B2*C3*ri3*t3;
calc=calc + 4*A4*si*t3;
calc=calc + 12*A3*B*si*t3;
calc=calc + A4*B*si*t3;
calc=calc + 12*A2*B2*si*t3;
calc=calc + 3*A3*B2*si*t3;
calc=calc + 4*A*B3*si*t3;
calc=calc + 3*A2*B3*si*t3;
calc=calc + A*B4*si*t3;
calc=calc + 16*A3*C*si*t3;
calc=calc + 36*A2*B*C*si*t3;
calc=calc + 5*A3*B*C*si*t3;
calc=calc + 24*A*B2*C*si*t3;
calc=calc + 12*A2*B2*C*si*t3;
calc=calc + 4*B3*C*si*t3;
calc=calc + 9*A*B3*C*si*t3;
calc=calc + 2*B4*C*si*t3;
calc=calc + 24*A2*C2*si*t3;
calc=calc + 36*A*B*C2*si*t3;
calc=calc + 9*A2*B*C2*si*t3;
calc=calc + 12*B2*C2*si*t3;
calc=calc + 15*A*B2*C2*si*t3;
calc=calc + 6*B3*C2*si*t3;
calc=calc + 16*A*C3*si*t3;
calc=calc + 12*B*C3*si*t3;
calc=calc + 7*A*B*C3*si*t3;
calc=calc + 6*B2*C3*si*t3;
calc=calc + 4*C4*si*t3;
calc=calc + 2*B*C4*si*t3;
calc=calc - 8*A4*ri*si*t3;
calc=calc - 22*A3*B*ri*si*t3;
calc=calc - 2*A4*B*ri*si*t3;
calc=calc - 18*A2*B2*ri*si*t3;
calc=calc - 6*A3*B2*ri*si*t3;
calc=calc - 2*A*B3*ri*si*t3;
calc=calc - 6*A2*B3*ri*si*t3;
calc=calc + 2*B4*ri*si*t3;
calc=calc - 2*A*B4*ri*si*t3;
calc=calc - 22*A3*C*ri*si*t3;
calc=calc + 2*A4*C*ri*si*t3;
calc=calc - 44*A2*B*C*ri*si*t3;
calc=calc + 4*A3*B*C*ri*si*t3;
calc=calc - 22*A*B2*C*ri*si*t3;
calc=calc - 3*A2*B2*C*ri*si*t3;
calc=calc - 8*A*B3*C*ri*si*t3;
calc=calc - 3*B4*C*ri*si*t3;
calc=calc - 18*A2*C2*ri*si*t3;
calc=calc + 6*A3*C2*ri*si*t3;
calc=calc - 22*A*B*C2*ri*si*t3;
calc=calc + 11*A2*B*C2*ri*si*t3;
calc=calc - 4*B2*C2*ri*si*t3;
calc=calc - 3*B3*C2*ri*si*t3;
calc=calc - 2*A*C3*ri*si*t3;
calc=calc + 6*A2*C3*ri*si*t3;
calc=calc + 8*A*B*C3*ri*si*t3;
calc=calc + 3*B2*C3*ri*si*t3;
calc=calc + 2*C4*ri*si*t3;
calc=calc + 2*A*C4*ri*si*t3;
calc=calc + 3*B*C4*ri*si*t3;
calc=calc + 4*A4*ri2*si*t3;
calc=calc + 8*A3*B*ri2*si*t3;
calc=calc + A4*B*ri2*si*t3;
calc=calc + 3*A3*B2*ri2*si*t3;
calc=calc - 8*A*B3*ri2*si*t3;
calc=calc + 3*A2*B3*ri2*si*t3;
calc=calc - 4*B4*ri2*si*t3;
calc=calc + A*B4*ri2*si*t3;
calc=calc + 2*A3*C*ri2*si*t3;
calc=calc - 4*A4*C*ri2*si*t3;
calc=calc - 4*A2*B*C*ri2*si*t3;
calc=calc - 18*A3*B*C*ri2*si*t3;
calc=calc - 14*A*B2*C*ri2*si*t3;
calc=calc - 19*A2*B2*C*ri2*si*t3;
calc=calc - 8*B3*C*ri2*si*t3;
calc=calc - 3*A*B3*C*ri2*si*t3;
calc=calc + 2*B4*C*ri2*si*t3;
calc=calc - 8*A2*C2*ri2*si*t3;
calc=calc - 13*A3*C2*ri2*si*t3;
calc=calc - 18*A*B*C2*ri2*si*t3;
calc=calc - 25*A2*B*C2*ri2*si*t3;
calc=calc - 10*B2*C2*ri2*si*t3;
calc=calc - 13*A*B2*C2*ri2*si*t3;
calc=calc + 2*A2*B2*C2*ri2*si*t3;
calc=calc - B3*C2*ri2*si*t3;
calc=calc - 6*A*C3*ri2*si*t3;
calc=calc - 8*A2*C3*ri2*si*t3;
calc=calc - 6*B*C3*ri2*si*t3;
calc=calc - 8*A*B*C3*ri2*si*t3;
calc=calc - A2*B*C3*ri2*si*t3;
calc=calc - 2*B2*C3*ri2*si*t3;
calc=calc + A*C4*ri2*si*t3;
calc=calc + B*C4*ri2*si*t3;
calc=calc + 2*A3*B*ri3*si*t3;
calc=calc + 6*A2*B2*ri3*si*t3;
calc=calc + 6*A*B3*ri3*si*t3;
calc=calc + 2*B4*ri3*si*t3;
calc=calc + 4*A3*C*ri3*si*t3;
calc=calc + 2*A4*C*ri3*si*t3;
calc=calc + 12*A2*B*C*ri3*si*t3;
calc=calc + 9*A3*B*C*ri3*si*t3;
calc=calc + 12*A*B2*C*ri3*si*t3;
calc=calc + 10*A2*B2*C*ri3*si*t3;
calc=calc + 4*B3*C*ri3*si*t3;
calc=calc + 2*A*B3*C*ri3*si*t3;
calc=calc - B4*C*ri3*si*t3;
calc=calc + 2*A2*C2*ri3*si*t3;
calc=calc + 7*A3*C2*ri3*si*t3;
calc=calc + 4*A*B*C2*ri3*si*t3;
calc=calc + 5*A2*B*C2*ri3*si*t3;
calc=calc + 2*B2*C2*ri3*si*t3;
calc=calc - 2*A*B2*C2*ri3*si*t3;
calc=calc - 2*A2*B2*C2*ri3*si*t3;
calc=calc - 2*B3*C2*ri3*si*t3;
calc=calc - 2*A*C3*ri3*si*t3;
calc=calc - 4*A2*C3*ri3*si*t3;
calc=calc - 4*A*B*C3*ri3*si*t3;
calc=calc - 2*A2*B*C3*ri3*si*t3;
calc=calc - B2*C3*ri3*si*t3;
calc=calc - 2*A4*si2*t3;
calc=calc - 2*A3*B*si2*t3;
calc=calc - 2*A4*B*si2*t3;
calc=calc + 2*A2*B2*si2*t3;
calc=calc - 3*A3*B2*si2*t3;
calc=calc + 2*A*B3*si2*t3;
calc=calc - A2*B3*si2*t3;
calc=calc - 8*A3*C*si2*t3;
calc=calc - 6*A2*B*C*si2*t3;
calc=calc - 9*A3*B*C*si2*t3;
calc=calc + 4*A*B2*C*si2*t3;
calc=calc - 15*A2*B2*C*si2*t3;
calc=calc + 2*B3*C*si2*t3;
calc=calc - 7*A*B3*C*si2*t3;
calc=calc - 12*A2*C2*si2*t3;
calc=calc - 6*A*B*C2*si2*t3;
calc=calc - 15*A2*B*C2*si2*t3;
calc=calc + 2*B2*C2*si2*t3;
calc=calc - 18*A*B2*C2*si2*t3;
calc=calc - 3*B3*C2*si2*t3;
calc=calc - 8*A*C3*si2*t3;
calc=calc - 2*B*C3*si2*t3;
calc=calc - 11*A*B*C3*si2*t3;
calc=calc - 6*B2*C3*si2*t3;
calc=calc - 2*C4*si2*t3;
calc=calc - 3*B*C4*si2*t3;
calc=calc + 4*A4*ri*si2*t3;
calc=calc + 2*A3*B*ri*si2*t3;
calc=calc + 4*A4*B*ri*si2*t3;
calc=calc - 8*A2*B2*ri*si2*t3;
calc=calc + 5*A3*B2*ri*si2*t3;
calc=calc - 6*A*B3*ri*si2*t3;
calc=calc - A*B4*ri*si2*t3;
calc=calc + 8*A3*C*ri*si2*t3;
calc=calc - A4*C*ri*si2*t3;
calc=calc - 4*A2*B*C*ri*si2*t3;
calc=calc + 6*A3*B*C*ri*si2*t3;
calc=calc - 18*A*B2*C*ri*si2*t3;
calc=calc + 17*A2*B2*C*ri*si2*t3;
calc=calc - 6*B3*C*ri*si2*t3;
calc=calc + 8*A*B3*C*ri*si2*t3;
calc=calc + A2*B3*C*ri*si2*t3;
calc=calc - B4*C*ri*si2*t3;
calc=calc - 3*A3*C2*ri*si2*t3;
calc=calc - 14*A*B*C2*ri*si2*t3;
calc=calc + 7*A2*B*C2*ri*si2*t3;
calc=calc - 10*B2*C2*ri*si2*t3;
calc=calc + 13*A*B2*C2*ri*si2*t3;
calc=calc - 2*A2*B2*C2*ri*si2*t3;
calc=calc + 2*B3*C2*ri*si2*t3;
calc=calc - 8*A*C3*ri*si2*t3;
calc=calc - 3*A2*C3*ri*si2*t3;
calc=calc - 8*B*C3*ri*si2*t3;
calc=calc + 3*A*B*C3*ri*si2*t3;
calc=calc + B2*C3*ri*si2*t3;
calc=calc - 4*C4*ri*si2*t3;
calc=calc - A*C4*ri*si2*t3;
calc=calc - 2*B*C4*ri*si2*t3;
calc=calc - 2*A4*ri2*si2*t3;
calc=calc + 2*A3*B*ri2*si2*t3;
calc=calc - 2*A4*B*ri2*si2*t3;
calc=calc + 10*A2*B2*ri2*si2*t3;
calc=calc - A3*B2*ri2*si2*t3;
calc=calc + 6*A*B3*ri2*si2*t3;
calc=calc + 3*A2*B3*ri2*si2*t3;
calc=calc + 2*A*B4*ri2*si2*t3;
calc=calc + 2*A3*C*ri2*si2*t3;
calc=calc + 2*A4*C*ri2*si2*t3;
calc=calc + 12*A2*B*C*ri2*si2*t3;
calc=calc + 8*A3*B*C*ri2*si2*t3;
calc=calc + 14*A*B2*C*ri2*si2*t3;
calc=calc - 2*A2*B2*C*ri2*si2*t3;
calc=calc + 4*B3*C*ri2*si2*t3;
calc=calc - 2*A*B3*C*ri2*si2*t3;
calc=calc - 2*A2*B3*C*ri2*si2*t3;
calc=calc + B4*C*ri2*si2*t3;
calc=calc + 10*A2*C2*ri2*si2*t3;
calc=calc + 5*A3*C2*ri2*si2*t3;
calc=calc + 14*A*B*C2*ri2*si2*t3;
calc=calc - 2*A2*B*C2*ri2*si2*t3;
calc=calc + 8*B2*C2*ri2*si2*t3;
calc=calc + B3*C2*ri2*si2*t3;
calc=calc + 6*A*C3*ri2*si2*t3;
calc=calc + A2*C3*ri2*si2*t3;
calc=calc + 4*B*C3*ri2*si2*t3;
calc=calc + 2*A*B*C3*ri2*si2*t3;
calc=calc + 2*A2*B*C3*ri2*si2*t3;
calc=calc - B2*C3*ri2*si2*t3;
calc=calc - 2*A*C4*ri2*si2*t3;
calc=calc - B*C4*ri2*si2*t3;
calc=calc - 2*A3*B*ri3*si2*t3;
calc=calc - 4*A2*B2*ri3*si2*t3;
calc=calc - A3*B2*ri3*si2*t3;
calc=calc - 2*A*B3*ri3*si2*t3;
calc=calc - 2*A2*B3*ri3*si2*t3;
calc=calc - A*B4*ri3*si2*t3;
calc=calc - 2*A3*C*ri3*si2*t3;
calc=calc - A4*C*ri3*si2*t3;
calc=calc - 2*A2*B*C*ri3*si2*t3;
calc=calc - 5*A3*B*C*ri3*si2*t3;
calc=calc + A*B3*C*ri3*si2*t3;
calc=calc + A2*B3*C*ri3*si2*t3;
calc=calc + 2*A2*C2*ri3*si2*t3;
calc=calc - 2*A3*C2*ri3*si2*t3;
calc=calc + 6*A*B*C2*ri3*si2*t3;
calc=calc + 10*A2*B*C2*ri3*si2*t3;
calc=calc + 5*A*B2*C2*ri3*si2*t3;
calc=calc + 2*A2*B2*C2*ri3*si2*t3;
calc=calc + 4*A*C3*ri3*si2*t3;
calc=calc + 8*A2*C3*ri3*si2*t3;
calc=calc + 3*A*B*C3*ri3*si2*t3;
calc=calc + A2*B*C3*ri3*si2*t3;
calc=calc - 2*A3*B*si3*t3;
calc=calc + A4*B*si3*t3;
calc=calc - 2*A2*B2*si3*t3;
calc=calc - 6*A2*B*C*si3*t3;
calc=calc + 4*A3*B*C*si3*t3;
calc=calc - 4*A*B2*C*si3*t3;
calc=calc + 3*A2*B2*C*si3*t3;
calc=calc - A*B3*C*si3*t3;
calc=calc - A2*B3*C*si3*t3;
calc=calc - 6*A*B*C2*si3*t3;
calc=calc + 6*A2*B*C2*si3*t3;
calc=calc - 2*B2*C2*si3*t3;
calc=calc + 3*A*B2*C2*si3*t3;
calc=calc - B3*C2*si3*t3;
calc=calc - 2*B*C3*si3*t3;
calc=calc + 4*A*B*C3*si3*t3;
calc=calc + B*C4*si3*t3;
calc=calc + 4*A3*B*ri*si3*t3;
calc=calc - 2*A4*B*ri*si3*t3;
calc=calc + 2*A2*B2*ri*si3*t3;
calc=calc + A3*B2*ri*si3*t3;
calc=calc - 2*A*B3*ri*si3*t3;
calc=calc + 2*A3*C*ri*si3*t3;
calc=calc + 12*A2*B*C*ri*si3*t3;
calc=calc - 5*A3*B*C*ri*si3*t3;
calc=calc + 4*A*B2*C*ri*si3*t3;
calc=calc - 5*A2*B2*C*ri*si3*t3;
calc=calc + 4*A*B3*C*ri*si3*t3;
calc=calc + 2*A2*B3*C*ri*si3*t3;
calc=calc + 6*A2*C2*ri*si3*t3;
calc=calc + 12*A*B*C2*ri*si3*t3;
calc=calc - 6*A2*B*C2*ri*si3*t3;
calc=calc + 2*B2*C2*ri*si3*t3;
calc=calc + 2*A*B2*C2*ri*si3*t3;
calc=calc + 2*A2*B2*C2*ri*si3*t3;
calc=calc + B3*C2*ri*si3*t3;
calc=calc + 6*A*C3*ri*si3*t3;
calc=calc + 4*B*C3*ri*si3*t3;
calc=calc - 2*A*B*C3*ri*si3*t3;
calc=calc + 2*B2*C3*ri*si3*t3;
calc=calc + 2*C4*ri*si3*t3;
calc=calc + B*C4*ri*si3*t3;
calc=calc - 2*A3*B*ri2*si3*t3;
calc=calc + A4*B*ri2*si3*t3;
calc=calc + 2*A2*B2*ri2*si3*t3;
calc=calc - 2*A3*B2*ri2*si3*t3;
calc=calc + 4*A*B3*ri2*si3*t3;
calc=calc - 2*A3*C*ri2*si3*t3;
calc=calc - 2*A2*B*C*ri2*si3*t3;
calc=calc + A3*B*C*ri2*si3*t3;
calc=calc + 6*A*B2*C*ri2*si3*t3;
calc=calc + 6*A2*B2*C*ri2*si3*t3;
calc=calc - 3*A*B3*C*ri2*si3*t3;
calc=calc - A2*B3*C*ri2*si3*t3;
calc=calc - 4*A2*C2*ri2*si3*t3;
calc=calc + A3*C2*ri2*si3*t3;
calc=calc + 8*A2*B*C2*ri2*si3*t3;
calc=calc - 5*A*B2*C2*ri2*si3*t3;
calc=calc - 2*A2*B2*C2*ri2*si3*t3;
calc=calc - 2*A*C3*ri2*si3*t3;
calc=calc + 2*A2*C3*ri2*si3*t3;
calc=calc - A*B*C3*ri2*si3*t3;
calc=calc - A2*B*C3*ri2*si3*t3;
calc=calc + A*C4*ri2*si3*t3;
calc=calc - 2*A2*B2*ri3*si3*t3;
calc=calc + A3*B2*ri3*si3*t3;
calc=calc - 2*A*B3*ri3*si3*t3;
calc=calc - 4*A2*B*C*ri3*si3*t3;
calc=calc - 6*A*B2*C*ri3*si3*t3;
calc=calc - 4*A2*B2*C*ri3*si3*t3;
calc=calc - 2*A2*C2*ri3*si3*t3;
calc=calc - A3*C2*ri3*si3*t3;
calc=calc - 6*A*B*C2*ri3*si3*t3;
calc=calc - 8*A2*B*C2*ri3*si3*t3;
calc=calc - 2*A*C3*ri3*si3*t3;
calc=calc - 4*A2*C3*ri3*si3*t3;
calc=calc + A4*t4;
calc=calc + 4*A3*B*t4;
calc=calc + 6*A2*B2*t4;
calc=calc + 4*A*B3*t4;
calc=calc + B4*t4;
calc=calc + 4*A3*C*t4;
calc=calc + 12*A2*B*C*t4;
calc=calc + 12*A*B2*C*t4;
calc=calc + 4*B3*C*t4;
calc=calc + 6*A2*C2*t4;
calc=calc + 12*A*B*C2*t4;
calc=calc + 6*B2*C2*t4;
calc=calc + 4*A*C3*t4;
calc=calc + 4*B*C3*t4;
calc=calc + C4*t4;
calc=calc - 2*A4*ri*t4;
calc=calc - 8*A3*B*ri*t4;
calc=calc - 12*A2*B2*ri*t4;
calc=calc - 8*A*B3*ri*t4;
calc=calc - 2*B4*ri*t4;
calc=calc - 6*A3*C*ri*t4;
calc=calc - 18*A2*B*C*ri*t4;
calc=calc - 18*A*B2*C*ri*t4;
calc=calc - 6*B3*C*ri*t4;
calc=calc - 6*A2*C2*ri*t4;
calc=calc - 12*A*B*C2*ri*t4;
calc=calc - 6*B2*C2*ri*t4;
calc=calc - 2*A*C3*ri*t4;
calc=calc - 2*B*C3*ri*t4;
calc=calc + A4*ri2*t4;
calc=calc + 4*A3*B*ri2*t4;
calc=calc + 6*A2*B2*ri2*t4;
calc=calc + 4*A*B3*ri2*t4;
calc=calc + B4*ri2*t4;
calc=calc + A3*C*ri2*t4;
calc=calc + 3*A2*B*C*ri2*t4;
calc=calc + 3*A*B2*C*ri2*t4;
calc=calc + B3*C*ri2*t4;
calc=calc - A2*C2*ri2*t4;
calc=calc - A3*C2*ri2*t4;
calc=calc - 2*A*B*C2*ri2*t4;
calc=calc - A2*B*C2*ri2*t4;
calc=calc - B2*C2*ri2*t4;
calc=calc - A*C3*ri2*t4;
calc=calc - A2*C3*ri2*t4;
calc=calc - B*C3*ri2*t4;
calc=calc + A3*C*ri3*t4;
calc=calc + 3*A2*B*C*ri3*t4;
calc=calc + 3*A*B2*C*ri3*t4;
calc=calc + B3*C*ri3*t4;
calc=calc + A2*C2*ri3*t4;
calc=calc + A3*C2*ri3*t4;
calc=calc + 2*A*B*C2*ri3*t4;
calc=calc + A2*B*C2*ri3*t4;
calc=calc + B2*C2*ri3*t4;
calc=calc - 2*A4*si*t4;
calc=calc - 6*A3*B*si*t4;
calc=calc - 6*A2*B2*si*t4;
calc=calc - 2*A*B3*si*t4;
calc=calc - 8*A3*C*si*t4;
calc=calc - 18*A2*B*C*si*t4;
calc=calc - 12*A*B2*C*si*t4;
calc=calc - 2*B3*C*si*t4;
calc=calc - 12*A2*C2*si*t4;
calc=calc - 18*A*B*C2*si*t4;
calc=calc - 6*B2*C2*si*t4;
calc=calc - 8*A*C3*si*t4;
calc=calc - 6*B*C3*si*t4;
calc=calc - 2*C4*si*t4;
calc=calc + 4*A4*ri*si*t4;
calc=calc + 11*A3*B*ri*si*t4;
calc=calc + 9*A2*B2*ri*si*t4;
calc=calc + A*B3*ri*si*t4;
calc=calc - B4*ri*si*t4;
calc=calc + 11*A3*C*ri*si*t4;
calc=calc + 22*A2*B*C*ri*si*t4;
calc=calc - 2*A3*B*C*ri*si*t4;
calc=calc + 11*A*B2*C*ri*si*t4;
calc=calc - 2*A2*B2*C*ri*si*t4;
calc=calc + 9*A2*C2*ri*si*t4;
calc=calc + 11*A*B*C2*ri*si*t4;
calc=calc - 2*A2*B*C2*ri*si*t4;
calc=calc + 2*B2*C2*ri*si*t4;
calc=calc + A*C3*ri*si*t4;
calc=calc - C4*ri*si*t4;
calc=calc - 2*A4*ri2*si*t4;
calc=calc - 4*A3*B*ri2*si*t4;
calc=calc + 4*A*B3*ri2*si*t4;
calc=calc + 2*B4*ri2*si*t4;
calc=calc - A3*C*ri2*si*t4;
calc=calc + 2*A2*B*C*ri2*si*t4;
calc=calc + 3*A3*B*C*ri2*si*t4;
calc=calc + 7*A*B2*C*ri2*si*t4;
calc=calc + 3*A2*B2*C*ri2*si*t4;
calc=calc + 4*B3*C*ri2*si*t4;
calc=calc + 4*A2*C2*ri2*si*t4;
calc=calc + 2*A3*C2*ri2*si*t4;
calc=calc + 9*A*B*C2*ri2*si*t4;
calc=calc + 2*A2*B*C2*ri2*si*t4;
calc=calc + 5*B2*C2*ri2*si*t4;
calc=calc + 3*A*C3*ri2*si*t4;
calc=calc + 2*A2*C3*ri2*si*t4;
calc=calc + 3*B*C3*ri2*si*t4;
calc=calc - A3*B*ri3*si*t4;
calc=calc - 3*A2*B2*ri3*si*t4;
calc=calc - 3*A*B3*ri3*si*t4;
calc=calc - B4*ri3*si*t4;
calc=calc - 2*A3*C*ri3*si*t4;
calc=calc - 6*A2*B*C*ri3*si*t4;
calc=calc - A3*B*C*ri3*si*t4;
calc=calc - 6*A*B2*C*ri3*si*t4;
calc=calc - A2*B2*C*ri3*si*t4;
calc=calc - 2*B3*C*ri3*si*t4;
calc=calc - A2*C2*ri3*si*t4;
calc=calc - 2*A3*C2*ri3*si*t4;
calc=calc - 2*A*B*C2*ri3*si*t4;
calc=calc - B2*C2*ri3*si*t4;
calc=calc + A*C3*ri3*si*t4;
calc=calc + A2*C3*ri3*si*t4;
calc=calc + A4*si2*t4;
calc=calc + A3*B*si2*t4;
calc=calc - A2*B2*si2*t4;
calc=calc - A3*B2*si2*t4;
calc=calc - A*B3*si2*t4;
calc=calc - A2*B3*si2*t4;
calc=calc + 4*A3*C*si2*t4;
calc=calc + 3*A2*B*C*si2*t4;
calc=calc - 2*A*B2*C*si2*t4;
calc=calc - A2*B2*C*si2*t4;
calc=calc - B3*C*si2*t4;
calc=calc + 6*A2*C2*si2*t4;
calc=calc + 3*A*B*C2*si2*t4;
calc=calc - B2*C2*si2*t4;
calc=calc + 4*A*C3*si2*t4;
calc=calc + B*C3*si2*t4;
calc=calc + C4*si2*t4;
calc=calc - 2*A4*ri*si2*t4;
calc=calc - A3*B*ri*si2*t4;
calc=calc + 4*A2*B2*ri*si2*t4;
calc=calc + 2*A3*B2*ri*si2*t4;
calc=calc + 3*A*B3*ri*si2*t4;
calc=calc + 2*A2*B3*ri*si2*t4;
calc=calc - 4*A3*C*ri*si2*t4;
calc=calc + 2*A2*B*C*ri*si2*t4;
calc=calc + 3*A3*B*C*ri*si2*t4;
calc=calc + 9*A*B2*C*ri*si2*t4;
calc=calc + 2*A2*B2*C*ri*si2*t4;
calc=calc + 3*B3*C*ri*si2*t4;
calc=calc + 7*A*B*C2*ri*si2*t4;
calc=calc + 3*A2*B*C2*ri*si2*t4;
calc=calc + 5*B2*C2*ri*si2*t4;
calc=calc + 4*A*C3*ri*si2*t4;
calc=calc + 4*B*C3*ri*si2*t4;
calc=calc + 2*C4*ri*si2*t4;
calc=calc + A4*ri2*si2*t4;
calc=calc - A3*B*ri2*si2*t4;
calc=calc - 5*A2*B2*ri2*si2*t4;
calc=calc - A3*B2*ri2*si2*t4;
calc=calc - 3*A*B3*ri2*si2*t4;
calc=calc - A2*B3*ri2*si2*t4;
calc=calc - A3*C*ri2*si2*t4;
calc=calc - 6*A2*B*C*ri2*si2*t4;
calc=calc - 4*A3*B*C*ri2*si2*t4;
calc=calc - 7*A*B2*C*ri2*si2*t4;
calc=calc + A2*B2*C*ri2*si2*t4;
calc=calc - 2*B3*C*ri2*si2*t4;
calc=calc - 5*A2*C2*ri2*si2*t4;
calc=calc - A3*C2*ri2*si2*t4;
calc=calc - 7*A*B*C2*ri2*si2*t4;
calc=calc + A2*B*C2*ri2*si2*t4;
calc=calc - 4*B2*C2*ri2*si2*t4;
calc=calc - 3*A*C3*ri2*si2*t4;
calc=calc - A2*C3*ri2*si2*t4;
calc=calc - 2*B*C3*ri2*si2*t4;
calc=calc + A3*B*ri3*si2*t4;
calc=calc + 2*A2*B2*ri3*si2*t4;
calc=calc + A*B3*ri3*si2*t4;
calc=calc + A3*C*ri3*si2*t4;
calc=calc + A2*B*C*ri3*si2*t4;
calc=calc + A3*B*C*ri3*si2*t4;
calc=calc - 2*A2*B2*C*ri3*si2*t4;
calc=calc - A2*C2*ri3*si2*t4;
calc=calc + A3*C2*ri3*si2*t4;
calc=calc - 3*A*B*C2*ri3*si2*t4;
calc=calc - 4*A2*B*C2*ri3*si2*t4;
calc=calc - 2*A*C3*ri3*si2*t4;
calc=calc - 2*A2*C3*ri3*si2*t4;
calc=calc + A3*B*si3*t4;
calc=calc + A2*B2*si3*t4;
calc=calc + A3*B2*si3*t4;
calc=calc + 3*A2*B*C*si3*t4;
calc=calc + 2*A*B2*C*si3*t4;
calc=calc + A2*B2*C*si3*t4;
calc=calc + 3*A*B*C2*si3*t4;
calc=calc + B2*C2*si3*t4;
calc=calc + B*C3*si3*t4;
calc=calc - 2*A3*B*ri*si3*t4;
calc=calc - A2*B2*ri*si3*t4;
calc=calc - 2*A3*B2*ri*si3*t4;
calc=calc + A*B3*ri*si3*t4;
calc=calc + A2*B3*ri*si3*t4;
calc=calc - A3*C*ri*si3*t4;
calc=calc - 6*A2*B*C*ri*si3*t4;
calc=calc - A3*B*C*ri*si3*t4;
calc=calc - 2*A*B2*C*ri*si3*t4;
calc=calc - 3*A2*C2*ri*si3*t4;
calc=calc - 6*A*B*C2*ri*si3*t4;
calc=calc - A2*B*C2*ri*si3*t4;
calc=calc - B2*C2*ri*si3*t4;
calc=calc - 3*A*C3*ri*si3*t4;
calc=calc - 2*B*C3*ri*si3*t4;
calc=calc - C4*ri*si3*t4;
calc=calc + A3*B*ri2*si3*t4;
calc=calc - A2*B2*ri2*si3*t4;
calc=calc + A3*B2*ri2*si3*t4;
calc=calc - 2*A*B3*ri2*si3*t4;
calc=calc - 2*A2*B3*ri2*si3*t4;
calc=calc + A3*C*ri2*si3*t4;
calc=calc + A2*B*C*ri2*si3*t4;
calc=calc + A3*B*C*ri2*si3*t4;
calc=calc - 3*A*B2*C*ri2*si3*t4;
calc=calc - 4*A2*B2*C*ri2*si3*t4;
calc=calc + 2*A2*C2*ri2*si3*t4;
calc=calc - 2*A2*B*C2*ri2*si3*t4;
calc=calc + A*C3*ri2*si3*t4;
calc=calc + A2*B2*ri3*si3*t4;
calc=calc + A*B3*ri3*si3*t4;
calc=calc + A2*B3*ri3*si3*t4;
calc=calc + 2*A2*B*C*ri3*si3*t4;
calc=calc + 3*A*B2*C*ri3*si3*t4;
calc=calc + 3*A2*B2*C*ri3*si3*t4;
calc=calc + A2*C2*ri3*si3*t4;
calc=calc + 3*A*B*C2*ri3*si3*t4;
calc=calc + 3*A2*B*C2*ri3*si3*t4;
calc=calc + A*C3*ri3*si3*t4;
calc=calc + A2*C3*ri3*si3*t4; return calc}