#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<windows.h>
#define CONSTVALUE 1000000000

void CFILE();
void FEHL(int LINDEX,int N,double Y0[],double YN[]);
void FUNC(double T,double X[],int N);
void DOUT();
void PARM();
void PLO(int I);
void RG2(double XE,double ABSERR,int NUM);
void RUNGE();
void SOLV(double T,double Y[]);
void INIT();
void OFILE();
int ICHEK(double i, int j,int k);
int INDEX(char a[],char b[]);
double DSIGN(double a,double b);
double FU(int I,double T,double X[]);
FILE *fp1,*fp2;
int NS,ING;
int ITC;
int ND;
int NOT;
int NOUT;
double ERR[130];
double T0,T1,TI;
extern double H,X[130],T;
extern double DX[130];
double DE[40],DF[40];
extern double OP[100];
double PT[4];
extern double PA[999];
extern double XX[100];
char LABEL[100][20];
int ERRNUM[130];
int QUOTIENT;
int ODD;

main(){
  int i,j=0,k=0;
  int l;
  int conum,count;

  T0=0;T1=0;TI=0;
  H=0;T=0;
  for(i=0;i<=129;i++)X[i]=0.0;
  for(i=0;i<=129;i++)DX[i]=0.0;
  for(i=0;i<=39;i++)DE[i]=0.0;
  for(i=0;i<=39;i++)DF[i]=0.0;
  for(i=0;i<=99;i++)OP[i]=0.0;
  for(i=0;i<=3;i++)PT[i]=0.0;
  for(i=0;i<=299;i++)PT[i]=0.0;
  for(i=0;i<=99;i++)XX[i]=0.0;
  OFILE();
  PARM();
  INIT();

  DOUT();
  PLO(j);
  j++;

  conum = (int)((CONSTVALUE*(double)QUOTIENT+(double)ODD)/(double)NOT);
  count = 0;

  for(l=0;l<ODD;l++){
	  RUNGE();
	  if(conum == count){
		DOUT();
		PLO(j);
		j++;
		count = 0;
	  }
	  count++;
  }

  for(l=0;l<QUOTIENT;l++){
	  for(i=1;i<=CONSTVALUE;i++){
		  RUNGE();
		  if(conum == count){
			  DOUT();
			  PLO(j);
			  j++;
			  count = 0;
		  }
		  count++;
	  }	  
  }

  CFILE();
  return(0);
}

void RUNGE(){
  int i,N;
  double CS1,CS2,X1[130],K0[130],PHI[130];
  N=NS+ING;
  for(i=0;i<=N-1;i++){
    X1[i]=X[i];
    K0[i]=0.0;
    DX[i]=0.0;
    PHI[i]=0.0;
  }
  CS1=1.0/sqrt(2.0);
  CS2=sqrt(2.0);
  
  SOLV(T,X1);
  FUNC(T,X1,N);
  for(i=0;i<=N-1;i++){
    X1[i]=X[i]+0.5*H*DX[i];
    PHI[i]+=DX[i];
    K0[i]=DX[i];
  }
  
  T+=0.5*H;
  SOLV(T,X1);
  FUNC(T,X1,N);
  for(i=0;i<=N-1;i++){
    X1[i]=X[i]+(CS1-0.5)*H*K0[i]+(1.0-CS1)*H*DX[i];
    PHI[i]+=(2.0-CS2)*DX[i];
    K0[i]=DX[i];
  }
  
  SOLV(T,X1);
  FUNC(T,X1,N);
  for(i=0;i<=N-1;i++){
    X1[i]=X[i]-CS1*H*K0[i]+(1.0+CS1)*H*DX[i];
    PHI[i]+=(2.0+CS2)*DX[i];
  }
  T+=0.5*H;
  SOLV(T,X1);
  FUNC(T,X1,N);
  for(i=0;i<=N-1;i++){
    PHI[i]+=DX[i];
    X[i]+=PHI[i]*H/6.0;
  }

  SOLV(T,X);
}

void SOLV(double T,double Y[]){
/* SOLVING CONSTRAINT EQUATIONS
   IF DEVIATION RATIO IS LESS THAN RELATIVE ERROR(D0)
   THEN THE PROGRAM OBTAINS SOLUTION */

  int i,I0,I1,ICHK;
  double DE1[40],DF1[40],D0,D1,D2;
  if(ND!=0){
    I0=ND/2;
    D0=1.0e-08;
/* D0 : RELATIVE ERROR */
    while(1){
      ICHK=0;
/* ICHK : INDICATOR FOR TERMINATION
        =1 ONE MORE REITERATION
        =0 TEMINATION */
      for(i=0;i<=I0-1;i++){
	D1=fabs(DE[i])*D0;
	I1=2*i-1;
	DE1[i]=FU(I1,T,Y);
	D2=fabs(DE1[i]-DE[i]);
	if(D2>D1) ICHK=1;
	DE[i]=DE1[i];
	
	D1=fabs(DF[i])*D0;
	if(D2>D1) ICHK=1;
	DF[i]=DF1[i];
      }
      if(ICHK==0) break;
    }
  }
}

void PLO(int I){
  int j,k;
  char CDIV[2]={','};

  if(I==0){
	  fprintf(fp2,"TIME         ,");
	  printf("TIME          ");
	  for(j=0;j<=NOUT-1;j++){
		  if(j==NOUT-1){
			  for(k=0;k<(signed int)(strlen(LABEL[j+1]));k++) 
				  if(LABEL[j+1][k]==' ') LABEL[j+1][k]='\0';
				  fprintf(fp2," %s", LABEL[j+1]);
				  printf(" %s",LABEL[j+1]);
		  }else{
			  for(k=0;k<(signed int)(strlen(LABEL[j+1]));k++) 
				  if(LABEL[j+1][k]==' ') LABEL[j+1][k]='\0';
				  fprintf(fp2,"%s%s",LABEL[j+1],CDIV);
				  printf(" %s%s",LABEL[j+1],CDIV);
		  }
		  if(!((j+1)%4)){
			  printf("\n");
			  printf("              ");
		  }
	  }
	  fprintf(fp2,"\n");
	  fprintf(fp2,"%e%s",T,CDIV);
	  printf("\n");
	  printf("%e%s",T,CDIV);
	  for(j=0;j<=NOUT-1;j++){
		  if(j==NOUT-1){
			  if(OP[j]>=0){
				  fprintf(fp2," %e",OP[j]);
				  printf(" %e",OP[j]);
			  }else{
				  fprintf(fp2,"%e",OP[j]);
				  printf("%e",OP[j]);
			  }
		  }else{
			  if(OP[j]>=0){
				  fprintf(fp2," %e%s",OP[j],CDIV);
				  printf(" %e%s",OP[j],CDIV);
			  }else{
				  fprintf(fp2,"%e%s",OP[j],CDIV);
				  printf("%e%s",OP[j],CDIV);
			  }
		  }
		  if(!((j+1)%4)){
			  printf("\n");
			  printf("              ");
		  }
	  }
	  fprintf(fp2,"\n");
	  printf("\n");
  }else{
	  fprintf(fp2,"%e%s",T,CDIV);
	  printf("%e%s",T,CDIV);
	  for(j=0;j<=NOUT-1;j++){
		  if(j==NOUT-1){
			  if(OP[j]>=0){
				  fprintf(fp2," %e",OP[j]);
				  printf(" %e",OP[j]);
			  }else{
				  fprintf(fp2,"%e",OP[j]);
				  printf("%e",OP[j]);
			  }
		  }else{
			  if(OP[j]>=0){
				  fprintf(fp2," %e%s",OP[j],CDIV);
				  printf(" %e%s",OP[j],CDIV);
			  }else{
				  fprintf(fp2,"%e%s",OP[j],CDIV);
				  printf("%e%s",OP[j],CDIV);
			  }
		  }
		  if(!((j+1)%4)){
			  printf("\n");
			  printf("              ");
		  }
	  }
	  fprintf(fp2,"\n");
	  printf("\n");
  }//end if(I==0)
}

void PARM(){
  int i,NO,N,I0;
  double V;
  char NA[4],NAM[20],LINE[72],DUMMY[20];
  
  NOUT=0;
  NOT=0;
  while(1){
    i=0;
    while((LINE[i]=fgetc(fp1))!='\n') i++;
    LINE[i]='\0';
    for(i=0;i<=2;i++) NA[i]=LINE[i];
    NA[i]='\0';
    if(strcmp(NA,"END")){
      if(!strcmp(NA,"PA ")){
	    for(i=0;i<=2;i++) DUMMY[i]=LINE[i+3];
	    DUMMY[i]='\0';
	    NO=atoi(DUMMY);
		i=0;
		while((DUMMY[i]=LINE[i+9])!=' ')i++;
	    DUMMY[i]='\0';
	    V=atof(DUMMY);
	    PA[NO]=V;							
      }else{
	    if(!strcmp(NA,"SU ")){
	      for(i=0;i<=2;i++) DUMMY[i]=LINE[i+2];
	      DUMMY[i]='\0';
	      NO=atoi(DUMMY);
		  i=0;
		  while((NAM[i]=LINE[i+6])!=' ')i++;
	      NAM[i]='\0';
		}else{
	      if(!strcmp(NA,"LA ")){
	        for(i=0;i<=2;i++) DUMMY[i]=LINE[i+2];
	        DUMMY[i]='\0';
	        NO=atoi(DUMMY);
			i=0;
		    while((NAM[i]=LINE[i+7])!=' ')i++;
	        for(i=i;i<=12;i++)NAM[i]=' ';
			NAM[i]='\0';
	        NOUT++;
	        strcpy(LABEL[NO],NAM);
		  }else{
	        if(!strcmp(NA,"NS ")){
	          for(i=0;i<=4;i++) DUMMY[i]=LINE[i+5];
	          DUMMY[i]='\0';
	          NO=atoi(DUMMY);
	          NS=NO;
			}else{
	          if(!strcmp(NA,"IN ")){
		        for(i=0;i<=4;i++) DUMMY[i]=LINE[i+5];
		        DUMMY[i]='\0';
		        NO=atoi(DUMMY);
		        ING=NO;
		        N=NS+ING;
		        for(i=1;i<=N;i++) X[i-1]=0.0e-00;
			  }else{
		        if(!strcmp(NA,"ND ")){
		          for(i=0;i<=4;i++) DUMMY[i]=LINE[i+5];
		          DUMMY[i]='\0';
		          NO=atoi(DUMMY);
		          ND=NO;
		          if(ND>0){
		            I0=ND/2;
		            for(i=1;i<=I0;i++){
		            DE[i]=0.0e-00;
		            DF[i]=0.0e-00;
					}
				  }
				}else{
		          if(!strcmp(NA,"PT ")){
		            for(i=0;i<=2;i++) DUMMY[i]=LINE[i+2];
		            DUMMY[i]='\0';
		            NO=atoi(DUMMY);
					i=0;
		            while((DUMMY[i]=LINE[i+8])!=' ')i++;
		            DUMMY[i]='\0';
		            V=atof(DUMMY);
		            PT[NO]=V;
				  }else{
		            if(!strcmp(NA,"NO ")){//Read Number of output times
		              for(i=0;i<11;i++) DUMMY[i]=LINE[i+5];
		              DUMMY[i]='\0';
		              NO=atoi(DUMMY);
		              NOT=NO;
					}else{
		              if(!strcmp(NA,"OP ")){
			            for(i=0;i<=11;i++) DUMMY[i]=LINE[i+5];
			            DUMMY[i]='\0';
			            NO=atoi(DUMMY);
			            NOUT=NO;
					  }else{
			            if(!strcmp(NA,"ST ")){
			              for(i=0;i<=2;i++) DUMMY[i]=LINE[i+2];
			              DUMMY[i]='\0';
			              NO=atoi(DUMMY);
						  i=0;
		                  while((DUMMY[i]=LINE[i+8])!=' ')i++;
			              DUMMY[i]='\0';
			              V=atof(DUMMY);
			              X[NO-1]=V;
						}else{
			              printf("      *** PARAMETER ERROR ***\n");
			              printf("      *** ILLEAGAL INPUT ***\n");
			              printf("      ...%s\n",LINE);
			              exit(1);
						}
					  }
					}
				  }
				}
			  }
			}
		  }
		}
      }
    }else{
      T0=PT[1];
      T1=PT[2];
      TI=PT[3];
      T=(T1-T0)/TI;
	  QUOTIENT = (int)(T/CONSTVALUE);
	  ODD = (int)(T - CONSTVALUE*(double)QUOTIENT);
      if(T <= NOT) NOT=(int)(T);
      H=TI;
      T=T0;
      break;
    }
  }
}

void INIT(){
  SOLV(T,X);
}

void OFILE(){
  char DYCP[50],NACP[50];
  printf("PLEASE INPUT DATA FILE NAME SUCH AS A:/DAT/SIM.BGS **\n");
  strcpy(NACP, "temp");
  strcpy(DYCP, NACP);

  if((fp1=fopen(strcat(DYCP,".PAR"),"rt"))==NULL){
    printf("FILE NOT OPEN\n");
    exit(1);
  }

  strcpy(DYCP,NACP);
  if((fp2=fopen(strcat(DYCP,".csv"),"wt"))==NULL){
    printf("FILE NOT OPEN\n");
    exit(1);
  }
}


void CFILE(){
  fclose(fp1);
  fclose(fp2);
}