//create the function to be integrated, sw is c, and sh is d

f := (1/(sw*sqrt(2*Pi)))*exp((-1*(x^2/(2*sw^2))))*(1/(sh*sqrt(2*Pi)))*exp((-1*(y^2/(2*sh^2))));
 

//Read in the data
line := readline(`means.txt`):
while line <> 0 do
    temp :=sscanf(",`%f %f %f %d`):
    w := temp[2]:
    ht := temp[3]:
    dist := temp[4]:
    ang := temp[1]:
 

    //use the appropriate c and d values, corresponding to the movement angle
    if (ang < 0.01) then g := subs(sw=evalf(.07165*dist),sh=evalf(.02836*dist),f) fi:
    if (ang < 0.4) and (ang > 0.3) then g := subs(sw=evalf(.06864*dist),sh=evalf(.03040*dist),f) fi:
    if (ang < 0.8) and (ang > 0.7) then g := subs(sw=evalf(.06342*dist),sh=evalf(.03439*dist),f) fi:
    if (ang < 1.2) and (ang > 1.1) then g := subs(sw=evalf(.05816*dist),sh=evalf(.03305*dist),f) fi:
    if (ang > 1.2) then g := subs(sw=evalf(.06654*dist),sh=evalf(.03454*dist),f) fi:
 

    //rotate the function by the movement angle (instead of rotating the target region)
    h := subs(x=xp*cos(ang)-yp*sin(ang),y=xp*sin(ang)+yp*cos(ang),g):

    //evaluate the integral
    j := evalf(int(int(h,yp=-ht/2..ht/2),xp=-w/2..w/2)):

    //print the results and get the next line of data
    printf(`%2.2f %2.2f %d %2.6f %2.10f\n`, w,ht,dist,ang,j);
    line := readline(`means.txt`):

od: