/*---------------------------------------------------------------------------*/
/* gccSURGE.c                                                                */
/*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define g 9.80665       /*重力加速度*/
#define NUMS 10         /*水槽断面設定最大点数*/
#define NUMQ 20         /*流量設定最大点数*/

char Title[256];        /*計算タイトル*/
int ICT;                /*計算ケ−ス設定*/
double AFCA;            /*AFC片側振幅流量(m3/s)*/
double AFCT;            /*AFC周期(s)*/

/*計算時間設定*/
double TMAX;            /*計算打切時間*/
double dt;              /*計算時間間隔*/
double DTWR;            /*打出時間間隔*/

/*制水口設定*/
double PAA;             /*制水口断面積*/
double PCI;             /*制水口流入時流量係数*/
double PCO;             /*制水口流出時流量係数*/

/*トンネル等設定*/
double RWL;             /*貯水池水位標高*/
double TNL;             /*トンネル延長*/
double TNA;             /*トンネル断面積*/
double TNC;             /*トンネル損失水頭係数*/

/*立坑断面*/
int NST;                /*入力断面数*/
double SAA[NUMS+2];     /*断面積*/
double SEL[NUMS+2];     /*下端標高*/
char SLB[NUMS+1][50];   /*立坑断面ラベル*/

/*流量設定*/
int NQT;                /*入力点数*/
double QTQ[NUMQ+1];     /*入力点流量*/
double QTI[NUMQ+1];     /*入力点時刻*/

/*その他*/
double AFCS;            /*AFC運転時の時刻指標*/
double QI;              /*初期流量*/

void DINP(char *fnameR);
void PRE();
int CALC(char *fnameW);
void WOUT1(FILE *fout);
void WOUT2(FILE *fout,double dts,double zw,double vw,double qn);
void RUNGE(double *qo,double vo,double zo,double fa,double dtt,double *dv,double *dz,double dts);
double FZ(double v,double q,double fa);
double FK(double v,double q,double cf);
double FV(double z,double v,double wrk);
double QNCAL(double tt);
/*---------------------------------------------------------------------------*/
 
int main(int argc,char *argv[])
{
    int i,nnn,nfile,iret;
    char *fnameR=argv[1]; /*入力ファイル名*/
    char *fnameW=argv[2]; /*出力ファイル名*/
    FILE *fp;

    DINP(fnameR);       /*データ入力*/
    PRE();              /*計算前処理*/
    iret=CALC(fnameW);  /*計算実行*/

    printf("%d : End of calculation.",iret);
    return 0;
}
/*---------------------------------------------------------------------------*/
void DINP(char *fnameR)
{
int i;
char dat[256];
char sv[20];
FILE *fin;
 
fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
        strncpy(Title,dat,strlen(dat)-1);    /*計算ケースタイトル*/
    fgets(dat,sizeof dat,fin);
        ICT=atoi(strtok(dat,","));           /*計算ケース設定(1:通常，2:AFC)*/
        AFCA=atof(strtok(NULL,","));         /*AFC片側振幅流量(m3/s)*/
        AFCT=atof(strtok(NULL,","));         /*AFC周期(s)*/
    fgets(dat,sizeof dat,fin);
        TMAX=atof(strtok(dat,","));          /*計算打切時間*/
        dt=atof(strtok(NULL,","));           /*計算時間間隔*/
        DTWR=atof(strtok(NULL,","));         /*打出時間間隔*/
    fgets(dat,sizeof dat,fin);
        PAA=atof(strtok(dat,","));           /*制水口断面積*/
        PCI=atof(strtok(NULL,","));          /*制水口流入時流量係数*/
        PCO=atof(strtok(NULL,","));          /*制水口流出時流量係数*/
    fgets(dat,sizeof dat,fin);
        RWL=atof(strtok(dat,","));           /*貯水池水位標高*/
        TNL=atof(strtok(NULL,","));          /*トンネル延長*/
        TNA=atof(strtok(NULL,","));          /*トンネル断面積*/
        TNC=atof(strtok(NULL,","));          /*トンネル損失水頭係数*/
    fgets(dat,sizeof dat,fin);
        NST=atoi(strtok(dat,","));           /*立坑入力断面数*/
        for(i=1;i<=NST;i++){
            fgets(dat,sizeof dat,fin);
            SAA[i]=atof(strtok(dat,","));    /*断面積*/
            SEL[i]=atof(strtok(NULL,","));   /*指定断面下端標高*/
            strcpy(sv,strtok(NULL,","));     /*立坑断面ラベル*/
            strncpy(SLB[i],sv,strlen(sv)-1); /*立坑断面ラベル*/
        }
    fgets(dat,sizeof dat,fin);
        NQT=atoi(strtok(dat,","));           /*流量設定点数*/
        for(i=1;i<=NQT;i++){
            fgets(dat,sizeof dat,fin);
            QTQ[i]=atof(strtok(dat,","));    /*入力点流量*/
            QTI[i]=atof(strtok(NULL,","));   /*入力点時刻*/
        }
fclose(fin);
 
}
/*---------------------------------------------------------------------------*/
void PRE()
{
    int i,j;
    double wa,wb;
    char wc[50];

    /*QI : 初期流量*/
    /*AFCS : 遮断開始時刻(通常計算では -1.0)*/
    /*通常計算ではQTI[1]に遮断開始時刻を入力*/
    /*AFC時計算ではQTI[2]に遮断開始時刻を入力*/

    QI=QTQ[1];
    switch(ICT){
        case 1:AFCS=-1.0  ;break;  /*通常*/
        case 2:AFCS=QTI[2];break;  /*AFC*/
    }

    /*立坑断面順番の確認(低標高順)*/
    for(i=1;i<=NST-1;i++){
        for(j=i+1;j<=NST;j++){
            if(SEL[i]>=SEL[j]){  /*標高の小さい順に並び替え*/
                wa=SEL[j];
                wb=SAA[j];
                strcpy(wc,SLB[j]);
                SEL[j]=SEL[i];
                SAA[j]=SAA[i];
                strcpy(SLB[j],SLB[i]);
                SEL[i]=wa;
                SAA[i]=wb;
                strcpy(SLB[i],wc);
            }
        }
    }
    SAA[NST+1]=SAA[NST];/*上端超過判定のためのダミー断面積*/
    SEL[NST+1]=SEL[NST];/*上端超過判定のためのダミー標高*/

}
/*---------------------------------------------------------------------------*/
int CALC(char *fnameW)
{
    int j,ic,np,npe,npn,npw,iw;
    double dto,qn,vi,zi,zo,tmi,tmx,zmi,zmx,dts,fa,dv,dz,zw,vw,dti,dtw;
    FILE *fout;

    fout=fopen(fnameW,"w");

    /*入力条件出力*/
    WOUT1(fout);

    np=(int)(DTWR/dt+0.00001); /*打出時間間隔内の計算step数*/
    dto=dt;                    /*計算時間間隔*/
    qn=QI;                     /*初期流量*/
    vi=qn/TNA;                 /*初期トンネル流速*/
    zi=TNC*vi*vi;              /*初期損失水頭*/
    if(vi<0.0)zi=-TNC*vi*vi;   /*初期損失水頭符号確定*/
    zo=RWL-zi;                 /*水槽内初期水位*/

    /*初期値出力（経過時間:0）*/
    dts=0.0;zw=zo;vw=vi;zi=zi;qn=qn;
    WOUT2(fout,dts,zw,vw,qn);

    /*計算開始断面検索*/
    for(j=NST;j>=1;j--){
        if(SEL[j]<zo)break;
    }

    npe=j;    /*計算開始断面番号*/
    npw=0;    /*打出時間間隔カウンタ初期化*/
    dts=0.0;  /*経過時間初期化*/

    do{
        dts=dts+dto;/*解を求めたい時刻*/
        npw=npw+1;
        fa=SAA[npe];
        RUNGE(&qn,vi,zi,fa,dto,&dv,&dz,dts);
        zw=RWL-zi-dz;

        iw=npe+1;
        if(iw>NST)iw=NST;
        if(zw<SEL[npe]||SEL[iw]<zw){
            dti=0.1*dto;
            for(ic=1;ic<=10;ic++){
                dtw=dts-dto+dti*(double)ic;
                zw=RWL-zi;
                if(SEL[npe+1]<zw)npe=npe+1;
                if(zw<SEL[npe])npe=npe-1;
                if(NST<npe){
                    /*USWL超過のため計算打ち切り*/
                    printf("***** Warning : USWL is NG. *****\n");
                    fprintf(fout,"#***** USWL is NG. *****\n");
                    fclose(fout);
                    return 1;
                }
                if(npe<1){
                    /*DSWL未満のため計算打ち切り*/
                    printf("***** Warning : DSWL is NG. *****\n");
                    fprintf(fout,"#***** DSWL is NG. *****\n");
                    fclose(fout);
                    return 2;
                }
                fa=SAA[npe];
                RUNGE(&qn,vi,zi,fa,dti,&dv,&dz,dtw);
                vi=vi+dv;
                zi=zi+dz;
            }
        }else{
            vi=vi+dv;
            zi=zi+dz;
        }
        /*打出時間間隔到達による出力*/
        if(npw>=np){
            npw=0;
            zw=RWL-zi;
            vw=vi;
            WOUT2(fout,dts,zw,vw,qn);
        }
     }while(dts<TMAX);

    fclose(fout);

    return 0;
}
/*---------------------------------------------------------------------------*/
void RUNGE(double *qo,double vo,double zo,double fa,double dtt,double *dv,double *dz,double dts)
{
    double wrk,dz0,dv0,cf,dtw,v1,z1,dv1,dz1,dv2,dz2,dv3,dz3;

    wrk=0.0;
    dz0=dtt*FZ(vo,*qo,fa);
    cf=PCI*PAA;              /*水位上昇(制水口流入)*/
    if(dz0>0.0)cf=PCO*PAA;   /*水位降下(制水口流出)*/
    wrk=FK(vo,*qo,cf);
    dv0=dtt*FV(zo,vo,wrk);

    v1=vo+0.5*dv0;
    z1=zo+0.5*dz0;
    dtw=dts-dtt*0.5;
    *qo=QNCAL(dtw);
    wrk=FK(v1,*qo,cf);
    dv1=dtt*FV(z1,v1,wrk);
    dz1=dtt*FZ(v1,*qo,fa);
 
    v1=vo+0.5*dv1;
    z1=zo+0.5*dz1;
    wrk=FK(v1,*qo,cf);
    dv2=dtt*FV(z1,v1,wrk);
    dz2=dtt*FZ(v1,*qo,fa);
 
    v1=vo+dv2;
    z1=zo+dz2;
    dtw=dts;
    *qo=QNCAL(dtw);
    wrk=FK(v1,*qo,cf);
    dv3=dtt*FV(z1,v1,wrk);
    dz3=dtt*FZ(v1,*qo,fa);
 
    *dv=(dv0+2.0*dv1+2.0*dv2+dv3)/6.0;
    *dz=(dz0+2.0*dz1+2.0*dz2+dz3)/6.0;
 
}
/*---------------------------------------------------------------------------*/
double FZ(double v,double q,double fa)
{
    /* Calculation of dz */
    return (q-TNA*v)/fa;
}
/*---------------------------------------------------------------------------*/
double FK(double v,double q,double cf)
{
    /* Calculation of k */
    return fabs((TNA*v-q)/cf)*(TNA*v-q)/cf/2.0/g;
}
/*---------------------------------------------------------------------------*/
double FV(double z,double v,double wrk)
{
    /* Calculation of dv */
    return (z-TNC*fabs(v)*v-wrk)/TNL*g;
}
/*---------------------------------------------------------------------------*/
double QNCAL(double tt)
{
    /* Calculation of discharge by interpolation method */
    int i;
    double qn;

    if(tt>=AFCS){
        for(i=2;i<=NQT;i++){
            if(QTI[i-1]<tt&&tt<=QTI[i]){
                qn=QTQ[i-1]+(QTQ[i]-QTQ[i-1])/(QTI[i]-QTI[i-1])*(tt-QTI[i-1]);
                break;
            }
        }
    }else{
        qn=QI+AFCA*sin(2.0*M_PI/AFCT*tt);
    }

    return(qn);
}
/*---------------------------------------------------------------------------*/
void WOUT1(FILE *fout)
{
    char dat[256];
    char sv[50];
    int i;
 
    strcpy(dat,"#,Title,");
    strcat(dat,Title);
    fprintf(fout,"%s\n",dat);
 
    fprintf(fout,"#\n");
    fprintf(fout,"#***********************\n");
    fprintf(fout,"# Input data \n");
    fprintf(fout,"#***********************\n");
 
    strcpy(dat,"#,ICT,AFCA,AFCT");
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#");
    sprintf(sv,"%d",ICT) ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",AFCA);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",AFCT);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
 
    strcpy(dat,"#,TMAX,dt,DTWR");
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#");
    sprintf(sv,"%.3f",TMAX);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",dt) ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",DTWR);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
 
    strcpy(dat,"#,PAA,PCI,PCO");
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#");
    sprintf(sv,"%.3f",PAA);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",PCI);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",PCO);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
 
    strcpy(dat,"#,RWL,TNL,TNA,TNC");
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#");
    sprintf(sv,"%.3f",RWL);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",TNL);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",TNA);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",TNC);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
 
    strcpy(dat,"#,NST");
    sprintf(sv,"%d",NST);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#,SHA,SHP,LBL");
    fprintf(fout,"%s\n",dat);
    for(i=1;i<=NST;i++){
        strcpy(dat,"#");
        sprintf(sv,"%.3f",SAA[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",SEL[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%s"  ,SLB[i]);strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
 
    strcpy(dat,"#,NQT");
    sprintf(sv,"%d",NQT);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#,QTE,TIM");
    fprintf(fout,"%s\n",dat);
    for(i=1;i<=NQT;i++){
        strcpy(dat,"#");
        sprintf(sv,"%.3f",QTQ[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",QTI[i]);strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
 
    fprintf(fout,"#\n");
    fprintf(fout,"#***********************\n");
    fprintf(fout,"# Result of calculation \n");
    fprintf(fout,"#***********************\n");
    strcpy(dat,"#Time,WL of Surge tank,Velocity of Tunnel,Discharge,k");
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"#(s),EL(m),(m/s),(m3/s),(m)");
    fprintf(fout,"%s\n",dat);
}
/*---------------------------------------------------------------------------*/
void WOUT2(FILE *fout,double dts,double zw,double vw,double qn)
{
    char dat[256];
    char sv[50];
    double k,cf;

    cf=0.5*(PCI+PCO)*PAA;
    k=FK(vw,qn,cf);

    sprintf(sv,"%e",dts);strcpy(dat,sv);
    sprintf(sv,"%e",zw);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%e",vw);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%e",qn);strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%e",k);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
}
/*---------------------------------------------------------------------------*/