# 2D Numerical Integration Package Version 0.8

num_int_2d := proc(u_range,v_range,fcn)
	local i,j,var1,a,b,ran1,ran2,var2,v_lower,v_upper,num_bars,v_lower_fcn,v_upper_fcn,v_lower_val,v_upper_val,nu,nv,u_val,int_val,v_val,total,u_total,v_expn:
	num_bars := 10:
	if nargs = 4 then num_bars := args[4] fi:

	var1 := op(1,u_range):
	var2 := op(1,v_range):

	ran1 := op(2,u_range):
	ran2 := op(2,v_range):

	a := evalf(op(1,ran1)):
	b := evalf(op(2,ran1)):

	v_lower := op(1,ran2):
	v_upper := op(2,ran2):
#print(u_range,v_range):
#print(var1,var2):
#print(ran1,ran2):
#print(a,b):
#print(v_lower,v_upper):
	v_lower_fcn := unapply(v_lower,var1):
	v_upper_fcn := unapply(v_upper,var1):

	nu := num_bars-1:
	nv := num_bars-1:

#	u_val := array(1..nu+1):
#	int_val := array(1..nu+1):

	total := 0:
	for i from 0 to nu do
#	    printf(`i=%d:\n`,i):
	    u_val := evalf((i/nu)*a+(1 - (i/nu))*b):
#	    print(`u_val,i,a,b`,u_val,i,a,b):


	    v_lower_val := v_lower_fcn(u_val):
	    v_upper_val := v_upper_fcn(u_val):
	    v_expn := subs(var1=u_val,fcn):
	    u_total := 0.5*evalf(subs(var2=v_lower_val,v_expn)):
	    j:= 'j':
#	    print(`nv,v_lower,v_upper`,nv,v_lower_val,v_upper_val):
	    u_total := u_total + sum(subs(var2=evalf((j/nv)*v_lower_val+
			  		  (1 - j/nv)*v_upper_val),
			        v_expn),
			   j=1..nv-1):
	    u_total := u_total +0.5* evalf(subs(var2=v_upper_val,v_expn)):
	    if (i=0) or (i=nu) then
		u_total := .5*u_total:
	    fi:
#	    print(u_total):
	    total := u_total*((v_upper_val-v_lower_val)/nv) + total:
	od:

	RETURN((total*(b-a))/nu):

end:

printf(`\t\tNumerical Double Integration _ Version 1.0\n \n`):
printf(`The function num_int_2d uses the trapezoidal rule in each\n`):
printf(`direction to numerically calculate double integrals.\n \n`):
printf(`Basic usage is, for example:\n \n`):
printf(`\tnum_int_2d(x = 1 .. 2, y = x .. x^2, sin(x+y));\n \n`):
printf(`to integrate sin(x+y) over the region 1 <=x<=2 and x<=y<=x^2\n`):
printf(`with a default of 10 divisions per integral.\n \n`):
printf(`Or:\n \n`):
printf(`\tnum_int_2d(x = 1 .. 2, y = x .. x^2, sin(x+y), 20);\n \n`):
printf(`to change the number of divisions per integral from 10 to 20.\n \n`):