/*
 * heat2d-e-sw-h.c -- solve the heat equation in a two dimensional square region
 * by explicit finite difference method. (version 2.0)
 *
 *	To compile this program  ->  ccmg heat2d-e-sw-h.c
 */

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

/* to use matrix, new_matrix() */
#include <matrix.h>

/* to use GLSC */
#define G_DOUBLE
#include <glsc.h>


int main()
{
    int             N, i, j, n, skip, nMax;
    double          d,b;
    double          h, lambda, tau, Tmax, t, c, dt ,yps1,yps2,ypsr,minyps,maxtau;
    double          x,r,hl,hr,hu,hd;
    double          xi,xil,xir,yj,yjd,yju;
    double          ul,ur,uu,ud; 
    matrix          u, newu;
    double          f(double, double, double, double);
    double          pita(double,double);
    double          a(double, double);

    printf("N: "); scanf("%d",&N);
    N=2*N;

    h = 1.0 / N;
    r=0.5;
    yps1=1.0;
    yps2=1.0;
    maxtau=1.0;
    minyps=1.0;
    ypsr=1.0e-10;

    d=0.3;
    b=1.05*d;

     /****************最小εと最大τを求める *****************/
      for (i = 1;i<N/2; i++){
	xi =r-i*h;
	xil=r-(i-1)*h;
	xir=r-(i+1)*h;

	for (j = 1;j<N/2; j++){
	  yj=r-j*h;
	  if(pita(xi,yj)<b*b-ypsr){
	    if(pita(xil,yj)>b*b-ypsr){

	      x=sqrt(d*d-yj*yj+sqrt(b*b*b*b-4*d*d*yj*yj));
	    
	      hl=x-xi;
	      yps1=hl/h;

	      if(yps1<minyps)
		minyps=yps1;
	      

	      x=sqrt(-d*d-xi*xi+sqrt(b*b*b*b+4*d*d*xi*xi));
	      hd=x-yj;

	      if(hd>h){
		hd=h;
	      }
	    
	    yps2=hd/h;

	    if(maxtau>yps1*yps2*h*h/(2.0*(yps1+yps2)))
	      maxtau=yps1*yps2*h*h/(2.0*(yps1+yps2));
	    }
	    if(pita(xir,yj)>b*b-ypsr){

	      x=sqrt(d*d-yj*yj-sqrt(b*b*b*b-4*d*d*yj*yj));
	    
	      hr=xi-x;
	      yps1=hr/h;

	      if(yps1<minyps)
		minyps=yps1;
	      

	      x=sqrt(-d*d-xi*xi+sqrt(b*b*b*b+4*d*d*xi*xi));
	      hd=x-yj;

	      if(hd>h){
		hd=h;
	      }
	    
	    yps2=hd/h;

	    if(maxtau>yps1*yps2*h*h/(2.0*(yps1+yps2)))
	      maxtau=yps1*yps2*h*h/(2.0*(yps1+yps2));
	    }
	  }
	}
      }
      printf("εの最小値=%g\n",minyps);

    if ((u = new_matrix(N + 1, N + 1)) == NULL) {
        fprintf(stderr, "配列 u を確保できませんでした。");
        exit(1);
    }
    if ((newu = new_matrix(N + 1, N + 1)) == NULL) {
        fprintf(stderr, "配列 newu を確保できませんでした。");
        exit(1);
    }
    printf("Tmax: "); scanf("%lf", &Tmax);

    printf("τ(≦%g ): ", maxtau);
    scanf("%lf", &tau);


    lambda = 2.0*tau/(h*h);
    printf("λ=%g になりました。\n", lambda);

    printf("Δt: ");
    scanf("%lf",  &dt);
    skip = rint(dt / tau);
    if (skip == 0) {
      printf("Δtが小さすぎるので、Δt=τとします。\n");
      skip = 1;
    }
    dt = skip * tau;

    g_init("Meta", 250.0, 160.0);
    g_device(G_BOTH);
    g_def_scale(0, 0.0, 2.0, 0.0, 2.0, 30.0, 70.0, 100.0, 72.0);
    g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
    g_sel_scale(0);

    for (i = 0; i <= N; i++)
        for (j = 0; j <= N; j++)
            u[i][j] = f(i*h, j*h, d, b);

    g_cls();
    g_hidden2(1.0, 1.0, 0.4, -1.0, 1.0, 5.0, 25.0, 20.0, 20.0, 20.0,
        150.0, 100.0, u, N + 1, N + 1,
        1, G_SIDE_NONE, 2, 1);

    nMax = rint(Tmax / tau);

    for (n = 1; n <= nMax; n++) {
      /**********************左半円***********************/
      for (i = 1; i < N/2; i++){
	xi =r-i*h;                          /* u[i][j]のx座標 */
	xil=r-(i-1)*h;                      /* u[i][j]の左の点のx座標*/
	xir=r-(i+1)*h;

	/**********************円盤領域左下*********************/
	    for (j = 1; j < N/2; j++){
	      yj=r-j*h;                     /* u[i][j]のy座標 */
	      yjd=r-(j-1)*h;                /* u[i][j]の下のy座標 */

	      if(pita(xi,yj) < b*b-ypsr){
		if(pita(xil,yj) > b*b-ypsr){
		  x=sqrt(d*d-yj*yj+sqrt(b*b*b*b-4*d*d*yj*yj));
		  hl=x-xi;
		  ul=a(r-x,j*h);
		}
		else{
		  hl=h;
		  ul=u[i-1][j];
		}
	      }
	      if(pita(xi,yj) < b*b-ypsr){
		if(pita(xir,yj) > b*b-ypsr){
		  x=sqrt(d*d-yj*yj-sqrt(b*b*b*b-4*d*d*yj*yj));
		  hr=xi-x;
		  ur=a(r-x,j*h);
		}
		else{
		  hr=h;
		  ur=u[i+1][j];
		}
	    }
	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xi,yjd) > b*b-ypsr){
		x=sqrt(-d*d-xi*xi+sqrt(b*b*b*b+4*d*d*xi*xi));
		hd=x-yj;
		ud=a(i*h,r-x);
	      }
	      else{
		hd=h;
		ud=u[i][j-1];
	      }
	      hu=h;
	      uu=u[i][j+1];

      	      newu[i][j]=u[i][j]+tau*(2.0*((ur-u[i][j])/hr-(u[i][j]-ul)/hl)/(hr+hl)+
		 	              2.0*((uu-u[i][j])/hu-(u[i][j]-ud)/hd)/(hu+hd));
	    }
	    else{
	      newu[i][j] = 0.0;
	    }
	    }

	    /**********************円盤領域左上**********************/
	  for (j = N/2; j < N; j++){
	    yj=j*h-r;
	    yju=(j+1)*h-r;
	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xil,yj) > b*b-ypsr){
		x=sqrt(d*d-yj*yj+sqrt(b*b*b*b-4*d*d*yj*yj));
		hl=x-xi;
		ul=a(r-x,j*h);
	      }
	      else{
		hl=h;
		ul=u[i-1][j];
	      }
	    }

	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xir,yj) > b*b-ypsr){
		x=sqrt(d*d-yj*yj-sqrt(b*b*b*b-4*d*d*yj*yj));
		hr=x;
		ur=a(r-x,j*h);
	      }
	      else{
		hr=h;
		ur=u[i+1][j];
	      }
	    }
	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xi,yju) > b*b-ypsr){
		x=sqrt(-d*d-xi*xi+sqrt(b*b*b*b+4*d*d*xi*xi));
		hu=x-yj;
		uu=a(i*h,r+x);
	      }
	      else{
		hu=h;
		uu=u[i][j+1];
	      }
	      hd=h;
	      ud=u[i][j-1];

       	      newu[i][j]=u[i][j]+tau*(2.0*((ur-u[i][j])/hr-(u[i][j]-ul)/hl)/(hr+hl)+
		         	      2.0*((uu-u[i][j])/hu-(u[i][j]-ud)/hd)/(hu+hd));
	    }
	    else  newu[i][j]=0.0;
	  }
      }
      /**********************右半円***********************/
      for (i = N/2; i < N; i++){
	xi=i*h-r;
	xir=(i+1)*h-r;
	xil=(i-1)*h-r;

	/***********************円盤領域右下**************************/
	  for (j = 1; j < N/2; j++){
	    yj=r-j*h;
	    yjd=r-(j-1)*h;
	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xir,yj) > b*b-ypsr){
		x=sqrt(d*d-yj*yj+sqrt(b*b*b*b-4*d*d*yj*yj));
		hr=x-xi;
		ur=a(r+x,j*h);
	      }
	      else{
		hr=h;
		ur=u[i+1][j];
	      }
	    }
	    if(pita(xi,yj )< b*b-ypsr){
	      if(pita(xil,yj) > b*b-ypsr){
		x=sqrt(d*d-yj*yj-sqrt(b*b*b*b-4*d*d*yj*yj));
		hl=xi-x;
		ul=a(r+x,j*h);
	      }
	      else{
	      hl=h;
	      ul=u[i-1][j];
	      }
	    }
	    if(pita( xi,yj) < b*b-ypsr){
	      if(pita( xi, yjd ) > b*b-ypsr){
		x=sqrt(-d*d-xi*xi+sqrt(b*b*b*b+4*d*d*xi*xi));
		hd=x-yj;
		ud=a(i*h,r-x);
	      }
	      else{
		hd=h;
		ud=u[i][j-1];
	      }
	      hu=h;
	      uu=u[i][j+1];
	      
	      newu[i][j]=u[i][j]+tau*(2.0*((ur-u[i][j])/hr-(u[i][j]-ul)/hl)/(hr+hl)+
				      2.0*((uu-u[i][j])/hu-(u[i][j]-ud)/hd)/(hu+hd));
	    }
	    else  newu[i][j] = 0.0;
	  }
	  /************************円盤領域右上***************************/
	  for (j = N/2; j < N; j++){
	    yj=j*h-r;
	    yju=(j+1)*h-r;
	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xir,yj) > b*b-ypsr){
		x=sqrt(d*d-yj*yj+sqrt(b*b*b*b-4*d*d*yj*yj));

		hr=x-xi;
		ur=a(r+x,j*h);
	      }
	      else{
		hr=h;
		ur=u[i+1][j];
	      }
	    }

	      if(pita(xi,yj) < b*b-ypsr){
		if(pita(xil,yj) > b*b-ypsr){
		  x=sqrt(d*d-yj*yj-sqrt(b*b*b*b-4*d*d*yj*yj));
		  hl=xi-x;
		  ul=a(r+x,j*h);
		}
		else{
		  hl=h;
		  ul=u[i-1][j];
		}
	      }
	    
	    if(pita(xi,yj) < b*b-ypsr){
	      if(pita(xi,yju) > b*b-ypsr){
		x=sqrt(-d*d-xi*xi+sqrt(b*b*b*b+4*d*d*xi*xi));
		hu=x-yj;
		uu=a(i*h,x+r);
	      }
	      else{
		hu=h;
		uu=u[i][j+1];
	      }
	      hd=h;
	      ud=u[i][j-1];
	      
	      newu[i][j]=u[i][j]+tau*(2.0*((ur-u[i][j])/hr-(u[i][j]-ul)/hl)/(hr+hl)+
				      2.0*((uu-u[i][j])/hu-(u[i][j]-ud)/hd)/(hu+hd));
	    }
	    else  newu[i][j]=0.0;
	  }
      }
        for (i = 0; i <= N; i++)
            for (j = 0; j <= N; j++)
                u[i][j] = newu[i][j];
        if (n % skip == 0){
            g_cls();
            g_hidden2(1.0, 1.0, 0.4, -1.0, 1.0, 5.0, 25.0, 20.0, 20.0, 20.0,
                150.0, 100.0, u, N + 1, N + 1,
                1, G_SIDE_NONE, 2, 1);
        }
        t = tau * n;
    }
    /* マウスでクリックされるのを待つ */
    g_sleep(-1.0);
    /* ウィンドウを消す */
    g_term();
    return 0;
}

double f(double x, double y,double d,double b){
  double r=0.5;
  if(sqrt(((x-r)*(x-r)+(y-r)*(y-r)+d*d)*((x-r)*(x-r)+(y-r)*(y-r)+d*d)-4*d*d*(x-r)*(x-r)) < b*b-1.0e-14 )
    return 15*(b*b-sqrt(((x-r)*(x-r)+(y-r)*(y-r)+d*d)*((x-r)*(x-r)+(y-r)*(y-r)+d*d)-4*d*d*(x-r)*(x-r)));

  else  return 0;
}

double pita(double x,double y){
  double d=0.3;
  double b=1.05*d;
  return sqrt((x*x+y*y+d*d)*(x*x+y*y+d*d)-4*d*d*x*x);
}

double a(double x, double y){
  return 0;
}