/******************************************************************************
        ソーンスウエイト法による月蒸発散量の計算
                                                              A.KONDOH 1992
 ******************************************************************************/

#include        <stdio.h>
#include        <stdlib.h>
#include        <fcntl.h>
#include        <sys\types.h>
#include        <sys\stat.h>
#include        <io.h>
#include        <math.h>
#include        <string.h>
#include        <ctype.h>
#include        <conio.h>
#include        <dos.h>
#include        <process.h>

#define         NUM             24   /* 月平均気温が26.5ﾟC以上の場合用のﾃｰﾌﾞﾙのﾃﾞ-ﾀｾｯﾄ数 */
#define         PI              3.141592653

#define         SE              stderr
#define         FILE_ERR        -1
#define         CRLF            puts("")
#define         CLS             printf("\033[2J")
#define         PAUSE           do{}while(!isspace(getch())) /* ﾎﾟｰｽﾞ */

/* ファイルディスクリプタ */

FILE    *fi;
FILE    *fo;

float   lat;        /* 緯度 */

float   Et[13];     /* ｿｰﾝｽｳｴｲﾄ法による蒸発散量 */
float   Ta[13];     /* 月平均気温 */
float   factor[13]; /* 補正係数 */

static int      nmonth[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};
int             nday[13][32];

static float    Th[24]={26.5, 27.0, 27.5, 28.0, 28.5, 29.0
                       ,29.5, 30.0, 30.5, 31.0, 31.5, 32.0
                       ,32.5, 33.0, 33.5, 34.0, 34.5, 35.0
                       ,35.5, 36.0, 36.5, 37.0, 37.5, 38.0};
static float    Eh[24]={13.50, 13.95, 14.37, 14.78, 15.17, 15.54
                       ,15.89, 16.21, 16.52, 16.80, 17.07, 17.31
                       ,17.53, 17.72, 17.90, 18.05, 18.18, 18.29
                       ,18.37, 18.43, 18.47, 18.49, 18.50, 18.50};

/* 関数ﾌﾟﾛﾄﾀｲﾌﾟ宣言 */

void    Thornthwate();
void    Over265(int);
void    get_input_parameters();
void    correct_Et();
float   Sun(int,int);
float   sekii(int,int);
void    cal_nday(int);
void    output_results();

/*-----------------------------------------------------------------------------
        ｿｰﾝｽｳｴｲﾄ法ﾒｲﾝﾌﾟﾛｸﾞﾗﾑ
  -----------------------------------------------------------------------------*/

void    main(void)
{

/* タイトル */

CLS;    fprintf(SE,"====== Evapotranspiration by Thornthwate Method ======\n");

/* データ入力 */

        get_input_parameters();

/* 月蒸発散量の計算 */

        Thornthwate();

/* 計算結果出力 */

        output_results();

}/* end of main */

/*-----------------------------------------------------------------------------
        ｿｰﾝｽｳｴｲﾄ法
  -----------------------------------------------------------------------------*/

void    Thornthwate()
{
        int        mon;
        double  power_a;
        double  ti;

        ti=0.0;
        for(mon=1;mon<=12;mon++){
                ti += pow( Ta[mon]/5.0 , 1.514);
        }

        power_a=(492390.0+17920.0*ti-77.1*ti*ti+0.675*ti*ti*ti)*1.0E-06;

        for(mon=1;mon<=12;mon++){

             /* 月平均気温が氷点下の場合 */

                if(Ta[mon]<0.0){
                        Et[mon]=0.0;
                        continue;
                }

             /* 月平均気温が26.5ﾟC以上の場合 */

                if(Ta[mon]>=26.5){
                        Over265(mon);
                        continue;
                }

             /* 普通の月 */

                Et[mon]=1.6*pow( (10.0*Ta[mon]/ti),power_a );

        }/* next month */

/* 昼の長さの補正 */

        correct_Et();

}/* end of Thornthwate */

/*-----------------------------------------------------------------------------
        月平均気温が26.5℃を越える月の処理（線型補間)
  -----------------------------------------------------------------------------*/

void    Over265(mon)
int     mon;
{
        int     i;
        float   answer;

        if(Ta[mon]>=38.0){
                Et[mon]=18.50;
        }else{
                for(i=0;i<NUM-1;i++){
                        if(Ta[mon]>=Th[i] && Ta[mon]<Th[i+1]){
                                answer=Eh[i]
                                      +((Eh[i+1]-Eh[i])*(Ta[mon]-Th[i])/(Th[i+1]-Th[i]));
                                break;
                        }
                }
                Et[mon]=answer;
        }
}

/*-----------------------------------------------------------------------------
        入力データ
  -----------------------------------------------------------------------------*/

void    get_input_parameters()
{
        int     year=2001;
        int     mon;
        char    fname[64];

CRLF;   fprintf(SE,"Input file name  : ");
        scanf("%s",fname);

        if((fi=fopen(fname,"rt"))==NULL){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(0);
        }
/*
        fprintf(SE,"Output file name : ");
        scanf("%s",fname);

        if((fo=fopen(fname,"wt"))==NULL){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(0);
        }
*/
/* データ読み込み */

        fscanf(fi,"%f",&lat);

        for(mon=1;mon<=12;mon++){
                fscanf(fi,"%f",&Ta[mon]);
        }

/* 月の日数の計算 */

        cal_nday(year);

}/* end of get_input_parameters */

/*-----------------------------------------------------------------------------
        月の長さの補正
  -----------------------------------------------------------------------------*/

void    correct_Et()
{
        int     mon;
        float   sd;           /* 可照時間 */

        for(mon=1;mon<=12;mon++){

             /* 可照時間の計算(15日の値で月を代表) */

                sd=Sun(mon,15);

             /* 係数の計算 */

                factor[mon]=(sd/12.0)*(nmonth[mon]/30.0);

             /* 補正 */

                Et[mon]=factor[mon]*Et[mon];

        }/* next month */

}/* end of correct_Et */

/*-------------------------------------------------------------------
        function Sun() 可照時間の計算
  -------------------------------------------------------------------*/

float   Sun(mon,day)
int     mon;
int     day;
{
        float   t1,t2,t0,cos_t0,sin_t0;
        float   sek,ret;
        float   lat2;             /* 緯度(ﾗｼﾞｱﾝ) */

        sek=sekii(mon,day);       /* 各月15日の赤緯 */

        lat2=lat*PI/180.0;

        cos_t0=-tan( (double)lat2 )*tan( (double)sek );
        sin_t0=sqrt( (double)( 1.0-cos_t0*cos_t0 ) );
        t0=atan2( (double)sin_t0,(double)cos_t0 );

        t1=-3.81972*t0+12.0;     /* 日出 */
        t2= 3.81972*t0+12.0;     /* 日没 */

        ret=(t2-t1);

        return ret;              /* 可照時間 */
}

/*-------------------------------------------------------------------
              function sekii() 赤緯
  -------------------------------------------------------------------*/

float   sekii(month,day)
int     month,day;
{
        float x,ret;
        x=(float)nday[month][day];            /* １月０日からの日数 */
        ret=0.01698*(x-80.0);                 /* ﾙﾝﾃﾞの式 */
        return( 0.4093*sin( (double)ret ) );  /* ﾙﾝﾃﾞの式 */
}

/*-------------------------------------------------------------------
        function cal_nday() 1月0日からの日数
  -------------------------------------------------------------------*/

void    cal_nday(year)
int     year;
{
        int     i,j,dummy;

        if( (year%4)==0 ){                           /* 閏年の処理 */
                if((year%100==0) && (year%400!=0)){  /* 但、100で割れるが400で */
                        nmonth[2] = 28;              /* 割れない年は平年 */
                }else{
                        nmonth[2] = 29;
                }
        }

        dummy=0;
        for(i=1;i<=12;i++){
                for(j=1;j<=nmonth[i];j++){
                        dummy++;
                        nday[i][j]=dummy;
                }
        }

}/* end of cal_nday */

/*-----------------------------------------------------------------------------
        計算結果の出力
  -----------------------------------------------------------------------------*/

void    output_results()
{
        int     mon;
        float   Total;

        Total=0.0;
CRLF;   for(mon=1;mon<=12;mon++){
                fprintf(SE,"Month : %2d  Ta : %6.1f  Et : %6.1f (factor:%4.2f)\n"
                          ,mon,Ta[mon],Et[mon]*10.0,factor[mon]);
                Total+=Et[mon]*10.0;
        }
        fprintf(SE,"Annual Et : %6.1f (mm)\n",Total);
}