#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "pssub.h"

double f(double x)
{
	 double fx;
     fx=sqrt(1.0-x*x);
     return fx;
}
double trapezoid(int n)
{
     int i; double dx,x,ans;
     x=-1.0; dx=2.0/(double)n;
     ans=f(x);
     for(i=1;i<=n-1;i++) {
       x=(double)i*dx-1.0;
       ans=ans+2.0*f(x);
      }
     x=1.0;
     ans=(ans+f(x))*dx/2.0;
     return ans;
}
void plotd(double *x1,double *y1,int n) 
{    int i;
     for (i=0;i<n; i++) { 
        x1[i]=log10(x1[i]);y1[i]=log10(y1[i]*100000); 
        printf("%lf %lf\n",x1[i],y1[i]);}
     linewidth(1.5);linety(1);
     setgray(0.0);
     plot(x1[0], y1[0], 3); 
     for (i=1;i<n;i++) {
      plot(x1[i],y1[i],2); }
     stroke();
     circ(x1[0],y1[0],0.03); stroke();
     for (i=1;i<n;i++) {
           circ(x1[i],y1[i],0.03); stroke();}
}
void frame()
{
     int ix,nx,iy,ny,ixx,iyy; double x,y,dx,dy; 

       linewidth(1.0);
       rect(0.0, 0.0, 4.0, 4.0);
       nx=4;
       for(ix=0;ix<=nx-1;ix++) {
         setgray(0.3);
         x = (double)ix; dy=4.0;
         if(ix!=0) {plot(x,0.0,3);plot(x,dy,2);stroke(); }
         for(ixx=2;ixx<=9;ixx++) {
             setgray(0.8);x = (double) ix+log10((double) ixx); 
             plot(x,0.0,3);plot(x,dy,2);stroke();}}
       ny=4;
       for(iy=0;iy<=ny-1;iy++) {
         setgray(0.5);
         y = (double)iy; dx=4.0;
         if(iy!=0) {plot(0.0,y,3);plot(dx,y,2);stroke(); }
         for(iyy=2;iyy<=9;iyy++) {
             setgray(0.8);y = (double) iy+log10((double) iyy); 
             plot(0.0,y,3);plot(dx,y,2);stroke();}}
      setgray(0.0);
      textx(3.5, 0.0, 1, "n");
        textx(0.0, 0.02, 1, "1"); textx(1.0, 0.02, 2, "10");
        textx(2.0, 0.02, 3, "100"); textx(3.0, 0.02, 4, "1000");
        textx(4.0, 0.02, 5, "10000");
      texty(0.0, 3.5, 3, "eps");
        texty(0.0, 0.0, 7, "0.00001");  texty(0.0, 1.0, 6, "0.0001");
        texty(0.0, 2.0, 5, "0.001");texty(0.0, 3.0, 4, "0.01");
        texty(0.0, 4.0, 3, "0.1");
     return;
}
void main()
{    
     int j,n,m=8;
     double exact,aint,aerr,rerr,x[8],y[8];
     exact=atan(1.0)*4.0/2.0;
     n=4;
     for (j=1;j<=8;j++) { 
        n=n*2;
        aint=trapezoid(n);
        aerr=fabs(aint-exact);
        rerr=aerr/exact;
        x[j-1]=(double) n;
        y[j-1]=rerr;
        printf("n= %3d, integral=%lf, error=%lf\n",n,aint,aerr);
      }
     init();
     viewport(0.2, 0.2, 0.8, 0.8); 
     xyworld(0.0, 0.0, 4.0, 4.0); 
     frame();
     plotd(x,y,m); 
     fin();
}