#! /usr/bin/calc

Fronteira=0; Disco = 0;
pi = 4*atan(1);
## Definição de P(x,y), Q(x,y), define DQ_Dx(x,y), define DP_Dy(x,y)
define P(x,y)		{return	re(-x*power(y,2));};
define Q(x,y)		{return re(-power(x,2)*y/2);};
define DQ_x(x,y)	{return re(-x*y);};
define DP_y(x,y)	{return re(-2*x*y);};
## define F(x,y) 		{return re(DQ_x(x,y)) -  re(DP_y(x,y));}

## as fronteiras do círculo  - fronteira do domínio da integral dupla
define h1(x,r,a,b) {return b + re(sqrt(power(r,2) - power(x-a,2))) ;} ;
define h2(x,r,a,b) {return b - re(sqrt(power(r,2) - power(x-a,2))) ;} ;
  

define riemann_disco(a,b, r, N){
	local soma=0, k=0;
	local Delta_x = 2*r/N, Delta_y=Delta_x;
	## local x=a-r,FimX=a+r, y, FimY;
	local x=a-r,y,FimY;
	while(k <= N-1){
		x = a-r + k*Delta_x;
		y = h2(x,r,a,b); FimY = h1(x,r,a,b); 
		if (y > FimY) swap(y,FimY);
		while(y < FimY){
			soma = soma + DQ_x(x,y) - DP_y(x,y)  ;
			y +=Delta_y;
		}
		k++;
	}
	return(soma*Delta_x*Delta_y);
}


## equações paramétricas da fronteira e suas derivadas 

## (x1(t), x2(t)) 
define x1(t,a,r){
	return re(a + r*cos(t));
}
define x2(t,b,r){              
        return re(b + r*sin(t));
}        

##  (dx1(t), x2(t))
define dx1(t, r){
        return re(-r*sin(t));
}
define dx2(t,r){
        return re(r*cos(t));       
}  

define riemann_fronteira(a,b, r, N){
	local soma=0,k=0;
	local t=0;
	local Delta_t = 2*pi/N;
	while( k <= N-1 ){
	t = k*Delta_t;
	soma = soma + P( x1(t,a,r),x2(t,b,r) )*dx1(t,r)  
		    + Q( x1(t,a,r),x2(t,b,r) )*dx2(t,r)  ;
	k++;
	}
	return(soma*Delta_t);
}

## Está errado!
## tipo matriz 2x1 de reais 
define 	riemann_paralela(a, b , r, N){
        local mat soma[] = {0,0}; ## disco, fronteira as integrais
	local Delta_t=2*pi/N;  			## intervalao de paramet [0, 2*pi] 
	local Delta_x=(2*r)/N,Delta_y=(2*r)/N;; ## integral dupla 
	local k=0, j;				## contadores das integrais
	local x, y, FimY;
	while(k <= N-1){
		soma[0] +=  P( x1(k*Delta_t,a,r),x2(k*Delta_t,b,r) )*dx1(k*Delta_t,r)
				  + Q( x1(k*Delta_t,a,r),x2(k*Delta_t,b,r) )*dx2(k*Delta_t,r);
		x = a-r + k*Delta_x;
		y = h2(x,r,a,b); FimY = h1(x,r,a,b); ## as partes da fronteira
		if (y > FimY) swap(y,FimY);
		## laço interno roda entre as curvas da fronteira 
		while(y < FimY){
			soma[1] += ( DQ_x(x,y) - DP_y(x,y)  )  ;
                        y +=Delta_y;
		}
		k++; ## controle do somatório externo
	}
	soma[0] = soma[0]*Delta_t;  ## integral na fronteira
	soma[1] = soma[1]*(Delta_x*Delta_y); ## integral no disco
	return(soma); ## tem duas coordenadas
}





define main(a, b, r, N){
	## local Fronteira, Disco;
	local mat soma[] = {0,0}; ## disco, fronteira:  as integrais
	soma = riemann_paralela(a, b, r, N); 
	Fronteira = soma[0]; Disco = soma[1];
	## Fronteira 	= riemann_fronteira(a,b, r, N);
	## Disco		= riemann_disco(a,b, r, N);
	print "Na fronteira = ",Fronteira;
	print "No disco     = ",Disco;
}



r = 1;
## delta = 0.005;
N = 1000;
a = 2; b = 1;
main(a, b, r, N);
## system("joe exer10_03.calc");