/* ---------------------------------- */
/* save as c18a1.c                    */
/* ---------------------------------- */
#include "x_hfile.h"
#include       "fa.h"
/* ---------------------------------- */
int main(void)
{
double ay = -3.;
double by =  3.;
int    ny = 2*100;
int    nx = 2*100;

double m =  0.;
/* ---------------------------------- */
 clrscrn();
 printf(" Let S be the part of the graph of z = %s with z >= 0.  \n\n", feq);
 printf(" Let C be the trace of S on the x-y-plane.              \n\n");
 printf(" Verify Stokes's theorem for the vector field,          \n\n");
 printf("     F(x,y,z) = %si %sj %sk\n\n\n\n\n",Meq,Neq,Peq);

 printf(" Stoke's theorem.\n\n");
 printf("     /           //                    \n");
 printf("    |           ||                     \n");
 printf("    O F.T ds  = || (curl F).n  dS      \n");
 printf("    |           ||                     \n");
 printf("   / C         //                      \n");
 printf("               S                   \n\n\n");
 stop();

/* ---------------------------------- */
 clrscrn();
 printf("     //                    \n");
 printf("    ||                     \n");
 printf("    || (curl F).n dS =     \n");
 printf("    ||                     \n");
 printf("   //                      \n");
 printf("   S                   \n\n\n");
 printf("   with F = Mi + Nj + Pk     \n\n");
 printf(" (curl F) = [(P_y-N_z)i + (M_y-P_z)j + (N_X-M_Y)k]\n\n\n\n");
 printf("            (-f_xi-f_yj+k)           \n");
 printf("        n =  ------------            \n");
 printf("           [(f_x)^2+(f_y)^2+1]^1/2 \n\n\n\n");
 printf("       dS = [(f_x)^2+(f_y)^2+1]^1/2 dA             dA = dxdy\n\n");
 stop();

/* ---------------------------------- */
 clrscrn();
 printf(" With the Stokes's theorem you find :\n\n\n");
 printf(" F : (x,y,z)-> %si %sj %sk \n\n",Meq,Neq,Peq);
 printf(" f :   (x,y)-> %s     \n\n", feq);
 printf(" u :     (y)-> %s       \n", ueq);
 printf(" v :     (y)-> %s     \n\n", veq);
 printf(" a = %+.1f   b = %+.1f\n\n",ay,by);
 
  m = stokes_dxdy(M,N,P,
                             f,
                             H,
                             ay,by,ny,
                             u,  v,nx);

 printf("     / b   / v(y)\n");
 printf("    |     |      \n");
 printf("    |     |   (curl F).n    [(f_x)^2+(f_y)^2+1]^1/2 dx dy = %.3f\n",m);
 printf("    |     |      \n");
 printf("   /  a  /   u(y)\n\n\n");
 stop();
 
 return 0;
}
/* ---------------------------------- */
/* ---------------------------------- */