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

#define   KE  201

int main( int argc, char **argv ){
  int    n,      k,     kc = KE/2,    NSTEPS = 230;
  int myrank,isize;
  double  *ex, *hy;
  double  ex_low_m1 = 0.0,  ex_low_m2 = 0.0;
  double  ex_high_m1 = 0.0, ex_high_m2 = 0.0;
  double  T = 0.0,  t0 = 40.0,   spread = 12.0, pulse;
  double t_str,elapsed_time;
  FILE   *fp, *fopen();
  
  MPI_Status stat;
  int s,f;
  int ilf,irt;
  int domain;
  double  *anse,*ansh; 
  
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
  MPI_Comm_size(MPI_COMM_WORLD,&isize);
  t_str = MPI_Wtime();
  
  //domain = (KE%isize==0 ? KE/isize : KE/isize+1);
  domain = ((KE-2)%isize==0 ? (KE-2)/isize : (KE-2)/isize+1)+1;//区間
  if(myrank==0)printf( "domain: %d \n",domain);
  
  s = (domain-1)*myrank;//始点
  f = (domain-1)*(myrank+1);//終点
  // printf( "rank %d %d : %d \n",myrank,s,f);

  ex = malloc(sizeof(double)*(domain+1));
  hy = malloc(sizeof(double)*(domain+1));
  if(myrank==0){
    //anse = malloc(sizeof(double)*KE);
    //ansh = malloc(sizeof(double)*KE);
    anse = malloc(sizeof(double)*(domain+1)*isize);
    ansh = malloc(sizeof(double)*(domain+1)*isize);
  }
    for ( k = 0; k < domain+1 ; k++ ){
        ex[k] = 0.0;
        hy[k] = 0.0;
    }

    ilf = myrank-1;
    if(myrank==0)ilf = MPI_PROC_NULL;
    irt = myrank+1;
    if(myrank==isize-1) irt = MPI_PROC_NULL;

    for ( n = 1 ; n <= NSTEPS ; n++ ){
      T = (double)n;

      for ( k = 1 ; k <domain  ; k++ ){
	ex[k] = ex[k] + 0.5 * ( hy[k-1] - hy[k] );
        }
      if(s<=kc && f>kc){
	  pulse  = ( t0 - T ) / spread;
	  pulse  = -0.5 * pulse * pulse;
	  pulse  = exp( pulse );
	  ex[kc-s+1] = pulse;
	}
     
      MPI_Sendrecv(&ex[1], 1, MPI_DOUBLE, ilf, 11,
		   &ex[domain], 1, MPI_DOUBLE, irt, 11,
		 MPI_COMM_WORLD, &stat);
     
      if(myrank==0){
        ex[0]    = ex_low_m2;
        ex_low_m2  = ex_low_m1;
        ex_low_m1  = ex[1];
	}
	if(myrank==isize-1){
        ex[ domain ] = ex_high_m2;
        ex_high_m2 = ex_high_m1;
        ex_high_m1 = ex[ domain-1 ];
	}
        for ( k = 0 ; k < domain ; k++ ){
	  hy[k] = hy[k] + 0.5 * ( ex[k] - ex[k+1] );
        }

	MPI_Sendrecv(&hy[domain-1], 1, MPI_DOUBLE, irt,12,
		     &hy[0], 1, MPI_DOUBLE, ilf, 12,
		     MPI_COMM_WORLD, &stat);
    }

    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Gather(ex+1,domain-1,MPI_DOUBLE,anse,domain-1,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Gather(hy+1,domain-1,MPI_DOUBLE,ansh,domain-1,MPI_DOUBLE,0,MPI_COMM_WORLD);

      free(ex);
      free(hy);

      elapsed_time = MPI_Wtime() - t_str;
      // printf( "Elapsed time : %f (msec)\n",1.0e3*elapsed_time );
    
if(myrank==0){
  //課題２用
  // fp = fopen( "time.dat", "a" );
  // fprintf(fp, "%f\n",1.0e3*elapsed_time );
  //  fclose( fp );

    fp = fopen( "Ex230.dat", "w" );
    for ( k = 0 ; k < KE ; k++ )
        fprintf( fp, "  %6.2f \n", anse[k] );
    fclose( fp );

    fp = fopen( "Hy230.dat", "w" );
    for ( k = 0 ; k < KE ; k++ )
        fprintf( fp, "  %6.2f \n", ansh[k] );
    fclose( fp );
    free(anse);
    free(ansh);
    } 
 
    return MPI_Finalize();
}