/******************************************************
*       Program to demonstrate two dimensional        *
*               Mueller's method                      *
* --------------------------------------------------- *
* Reference: BASIC Scientific Subroutines, Vol. II    *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].   *
*                                                     *
*                  C++ version by J-P Moreau, Paris.  *
*                         (www.jpmoreau.fr)           *
* --------------------------------------------------- *
* Example: Find a real root of f(x,y)=(x+1)^5*(y-1)^5 *
*                                                     *
* Sample run:                                         *
*                                                     *
* Input the initial guesses and their bonds:          *
*    X0 = 0                                           *
*    Bond on X0 = 3                                   *
*    Y0 = 0                                           *
*    Bond on Y0 = 3                                   *
* Error criterion: 1e-6                               *
* Maximum number of iterations: 100                   *
*                                                     *
* The calculated zero is X= -0.999998  Y= 0.999998    *
* The associated W value is W = -0.000000             *
* The number of steps was: 43                         *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

// global variables
int    k,n;
double b1,b2,e,x,y,x0,yy0;

// shared by Utility() and Mueller2D()
double a1,b,c1,d1,e1,e2,e3,xl,xl1;


//**********************************************
//  Function subroutine                          
double W(double x, double y) {
  double temp;
  temp = (x+1)*(x+1)*(x+1)*(x+1)*(x+1);
  return (temp * (y-1)*(y-1)*(y-1)*(y-1)*(y-1));
}
//**********************************************

  // Utility function used by Mueller2D()
  void Utility() {
    a1=xl1*xl1*e1-d1*d1*e2+(xl1+d1)*e3;
    c1=xl1*(xl1*e1-d1*e2+e3);
    b=a1*a1-4.0*d1*c1*e3;
    // Test for complex root, meaning the parabola is inverted
    if (b<0) b=0;
    // Choose closest root
    if (a1<0) a1=a1-sqrt(b);
    if (a1>0) a1=a1+sqrt(b);
    // Guard against a divide by zero
    if (fabs(a1)+fabs(b)==0) a1=4.0*d1*e3;
    // Calculate a relative distance of next gues
    // and guard against a divide by zero
    if (a1==0) a1=1e-7;
    xl=-2.0*d1*e3/a1;
  }

/************************************************
*  Two Dimensional Mueller's method subroutine  *
* --------------------------------------------- *
* This routine iteratively seeks the root of a  *
* function W(x,y) by fitting a parabola to      *
* three points and calculating the nearest root *
* as described in Becket and Hurt, numerical    *
* calculations and algorithms.                  *
* INPUTS:                                       *
*   x0, y0 - The initial guess                  *
*   b1, b2 - A bound on the error in this guess *
*   e - The convergence criteria                *
*   n - The maximum number of iterations        *
* OUTPUTS:                                      *
*   x,y - The value of the root found           *
*   k - The number of iterations performed.     *
************************************************/
  void Mueller2D() {
  //Label: e100
  double x1,x2,x3,y1,y2,y3;
  // Set up the three evaluation points
  k=1;
e100: x3=x0; x1=x3-b1; x2=x3+b1;
  // Calculate Mueller parameters and guard against divide by zero
  if (x2-x1==0) x2=x2*1.0000001;
  if (x2-x1==0) x2=x2+0.0000001;
  xl1=(x3-x2)/(x2-x1);
  d1=(x3-x1)/(x2-x1);
  // Get values of function
  e1=W(x1,yy0);
  e2=W(x2,yy0);
  e3=W(x3,yy0);
  Utility();
  // Calculate new x estimate
  b1=xl*(x3-x2);
  x=x3+b1;
  // Test for convergence
  if (fabs(b1)+fabs(b2)<e) return;
  x0=x;
  // Repeat for the y direction
  y3=yy0; y1=y3-b2; y2=y3+b2;
  // Calculate Mueller parameters and guard against divide by zero
  if (y2-y1==0) y2=y2*1.0000001;
  if (y2-y1==0) y2=y2+0.0000001;
  xl1=(y3-y2)/(y2-y1);
  d1=(y3-y1)/(y2-y1);
  // Get values of function
  e1=W(x0,y1);
  e2=W(x0,y2);
  e3=W(x0,y3);
  Utility();
  // Calculate new y estimate
  b2=xl*(y3-y2);
  y=y3+b2;
  // Test for convergence
  if (fabs(b1)+fabs(b2)<e) return;
  // Test for number of iterations
  if (k>=n) return;
  yy0=y; k++;
  // Start another pass
  goto e100;
} //Mueller2D()


void main()  {
  printf("\n Input the initial guess:\n\n");
  printf("    X0 = "); scanf("%lf",&x0);
  printf("    Bond on X0 = "); scanf("%lf",&b1);
  printf("\n    Y0 = "); scanf("%lf",&yy0);
  printf("    Bond on Y0 = "); scanf("%lf",&b2);
  printf("\n Error criterion: "); scanf("%lf",&e);
  printf(" Maximum number of iterations: "); scanf("%d",&n);

  Mueller2D();  // Call 2D Mueller's routine

  printf("\n\n The calculated zero is (X,Y) = %f %f\n", x, y);
  printf(" The associated Y value is Y = %f\n", W(x,y));
  printf(" The number of steps was: %d\n\n", k);
}

// End of file Mueller2.cpp