Question: Double Integral of a Trig function

Hello Everyone,

I would like  to write a code to calculate the riemann sum approximation of the curve sin(sqrt(x^2+y^2)+1 

calculate the actual volume using integration

use the ranges x=[-2pi to 2pi]=y, and 20 subdivisions.

also display the curve and the parallelepiped approximations on the same plot.

 my code give me  an error..please someone can advise me which part of my code is not correct..

Thanks very much in advance!

> g := proc (x, y) options operator, arrow; sin(sqrt(x^2+y^2))+1 end proc;

 Z := g(x, y);

 X := Array(0 .. n); Y := Array(0 .. n); Z := Array(0 .. n);

c := -2*Pi; d := 2*Pi;

a := -2*Pi; b := 2*Pi;

n := 15;

s := 0;

X[0] := a; Y[0] := c; Z := g[X[0], Y[0]]; s := 0;

dX := (b-a)/N;dY := (d-c)/N;

    for i from 0 to N do      for j from o to N do     s:=s+Z[i]*dX*dY     X[i]:=X[i-1]+dX: Y[i]:=Y[i-1]+dY: Z[i]:=g([Xi],[Yi])  end do;

This the error message i am getting (Error, unterminated loop)

>