/*------------------------------------------------------------------------------
        GIMMSデータのトレンド解析
  ------------------------------------------------------------------------------*/

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

#define         NUMBER_OF_LINES		2091
#define         NUMBER_OF_PIXELS        4950
#define         NUMBER_OF_CHANNELS	25

signed short    NDVI[NUMBER_OF_CHANNELS][NUMBER_OF_PIXELS];
float           gain[NUMBER_OF_PIXELS];

void    main(argc,argv)
int     argc;
char    *argv[];
{
        int	fi,fo;
	int	i,j,ch;
        int     start_channel,end_channel;
        float   sum_x,sum_y,num;
        float   ave_x,ave_y,Sx,Sy,Sxy;

	if(argc!=5){
		printf("USAGE: TREND [in(BIL)] [out] [start year] [end year]\n");
		exit(-1);
	}

        if( (fi=open(argv[1],O_RDONLY | O_BINARY))==-1){
                       printf("Can't open %s.\n",argv[1]);
                       exit(-1);
        }

        if( (fo=open(argv[2],O_RDWR | O_TRUNC | O_CREAT | O_BINARY,
                        S_IREAD | S_IWRITE))==-1){
                       printf("Can't open %s.\n",argv[2]);
                       exit(-1);
        }

        start_channel=atoi(argv[3])-1982;
        end_channel=atoi(argv[4])-1982;
        if(end_channel>NUMBER_OF_CHANNELS){
                printf("End channel error.\n");
                exit(-1);
        }

	for(i=0;i<NUMBER_OF_LINES;i++){
                if( !(i%100) ) printf("*");

                for(ch=0;ch<NUMBER_OF_CHANNELS;ch++){
                        read(fi,NDVI[ch],NUMBER_OF_PIXELS*sizeof(short));
                }

                for(j=0;j<NUMBER_OF_PIXELS;j++){

                        if(NDVI[0][j]==-9999){
                                gain[j]=-99999.9;
                                continue; /* 領域外 */
                        }

                        sum_x=sum_y=num=0.0;
                        for(ch=start_channel;ch<=end_channel;ch++){
                                if(NDVI[ch][j]>0){
                                        num+=1.0;
                                        sum_x+=ch;
                                        sum_y+=(float)NDVI[ch][j]/1000.0;
                        }
                        }
			ave_x=sum_x/num;
			ave_y=sum_y/num;

			Sx=Sy=Sxy=0.0;
			for(ch=start_channel;ch<=end_channel;ch++){
      				Sx=Sx+((float)ch-ave_x)*((float)ch-ave_x);
        			Sy=Sy+((float)NDVI[ch][j]/1000.0-ave_y)*((float)NDVI[ch][j]/1000.0-ave_y);
	        		Sxy=Sxy+((float)ch-ave_x)*((float)NDVI[ch][j]/1000.0-ave_y);
                        }

			gain[j]=Sxy/Sx;

		}

		write(fo,gain,NUMBER_OF_PIXELS*sizeof(float));
	}

        printf("\n");

}/* end of main */