/*
 * newton2.c -- Newton 法で方程式 f(x)=0 を解く
 *  コンパイル gcc -o newton2 newton2.c -lm
 *  いずれも実行可能ファイルの名前は newton で、実行は ./newton
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

double r;

int main ()
{
  int i,j, maxitr = 100;
  double f(double), dfdx(double), x, dx, eps;
  double pai2;

  printf("r=");
  scanf("%lf",&r);

  printf("j=");
  scanf("%d",&j);

//pai2=pai/2
  pai2=2*atan(1.0);

  if(j<1){
    printf("Please Input j>=1 int");
    exit(0);
  }

  //auto pick out x0
  if(r<pai2){
    if(j==1){
      x=(r+pai2)/2;
    }
    if(j>1){
      x=2*pai2*(j-1);
    }
  }
  if(r>pai2){
    if(j==1){
      x=pai2+0.0001;
    }
    if(j>1){
      x=2*pai2*(j-1);
    }
  }
  if(r==pai2){
    if(j==1){
      x=pai2+0.001;
    }
    if(j>1){
      x=2*pai2*(j-1);
    }
  }
  if(r<=0){
    printf("Input r>0\n");
    exit(0);
  }
  //accuracy(seido)
  eps=1e-15;

  for (i = 0; i < maxitr; i++) {
    dx = - f(x) / dfdx(x);
    x += dx;
    printf("f(%20.15f)=%9.2e\n", x, f(x));
    if (fabs(dx) <= eps)
      break;
  }
  printf("lambda_%d=%20.15f\n",j,x*x);

  //gnuplot

  FILE *hoge;

  hoge=popen("gnuplot -persist","w");
  fprintf(hoge,"set xrange[0:1]\n");
  fprintf(hoge,"plot cos(%f*x)+%f*sin(%f*x)/%f \n",x,r,x,x);

  pclose(hoge);

  return 0;
}


double f(double x)
{
  return tan(x) - 2*r*x/(x*x-r*r);
}

/* 関数 f の導関数 (df/dx のつもりで名前をつけた) */
double dfdx(double x)
{
  return 1/(cos(x)*cos(x))+2*r*(2*x*x+r*r)/((x*x-r*r)*(x*x-r*r));
}