Mann Kendall 程序

能对气候水文等相关数据序列进行突变检验（fortran）和趋势分析(C)。

//This program calculates Mann-Kendall parameters for trend test.
//Condition: x is ordered from low to high.
//These parameters are: 
//    n-1  n
// S=sum （sum （sgn（yj-yi））; 
//   i=1 j=i+1
// n: Sample size;
//sigma = sqrt（n（n-1）（2n+5）/18）
// Zs =（S-1）/sigma        S>0
// Zs = 0                 S=0
// Zs = （S+1）/sigma       S<0
// Get the critical value of standard normal distribution with a probability of exceedance of alph/2.
// Get the Mann-Kendal slope 

#include 
#include 
#include 

#define ELEM_SWAP（ab） { register float t=（a）;（a）=（b）;（b）=t; }

double median（double *arr int n） ;
double mk_magnitude（float *x double *y int n）;


main （int argc char *argv[]）  {

 int nrows i j k s positive negative equal;
 float   *x temp tau sigma zs;
 double  *y slope; 
 float norm[6]={0.8421.2821.6451.9602.3262.576};
 float slevel[6]={0.600.800.900.950.980.99}level_out;

 FILE *fi *fo;

 if（argc !=3） {
  printf（“USAGE %s input output\n“ argv[0]）;
  exit（0）;
 }

 if（（fi = fopen（argv[1] “r“）） == NULL） {
  printf（“Cannot open input file.\n“）;
  exit（0）;
 }

 if（（fo = fopen（argv[2] “w“）） == NULL） {
  printf（“Cannot create output file.\n“）;
  exit（0）;
 }

 nrows = GetLineNumbers（fi）;
 printf（“Row number is %d\n“ nrows）;

 x= （float *） calloc（nrows sizeof（float））;
 y = （double *） calloc（nrows sizeof（double））;

 for（i=0; i   fscanf（fi “%f %lf“ &x[i] &y[i]）;
//   printf（“%f %lf\n“ x[i] y[i]）;
 }

 positive = negative = equal=0; 
 for（i=0; i   temp = y[i];
   for（j=i+1; j      if（temp == y[j]） equal++;
      if（temp < y[j]） positive++;
      if（temp > y[j]） negative++;
   }
  }

  s= positive + negative *（-1）;
  tau = s/ （nrows * （nrows-1） /2）;
  sigma = sqrt（nrows * （nrows-1） * （2*nrows+5）/18）;
  if（s>0） zs= （s-1）/sigma;
  if（s==0） zs=0;
  if（s<0） zs= （s+1）/sigma;

    level_out=-9;
    for（i=0;i<6;i++） {
      if（fabs（zs）>=norm[i]）  {
	level_out=slevel[i];
      } 
    }
   
//    if（level_out != -9）  {
        slope= mk_magnitude（xy nrows）;
        printf（“done slope\n“）;
//        }
//    else slope =0;

  printf（“zs=%f s=%d positive= %d negative=%d equal=%d level_out=%f slope=%lf\n“ zs s positive negative equal level_out slope）;
  fprintf（fo “Zs=%f S=%d N=%d Sigma=%f\n“ zs s nrows sigma）;
  fprintf（fo “Positive=%d Negative=%d Equal=%f\n“ positive negative equal）;
  fprintf（fo “level_out=%f slope=%lf\n“ level_out slope）;

 }
 
 int GetLineNumbers（FILE *fp） {
  int nlines=0;
  char ch;
  long start;
  start = ftell（fp）;
  while（（ch=fgetc（fp）） != EOF） if（ch==‘\n‘） nlines++;
  fseek（fp start SEEK_SET）;
  return nlines;
}

double mk_magnitude（float *x double *y int n）
{
  double *D;
  double b;
  int i j l;
  int nd;
  
  /* Calculate trend magnitude */
  nd = （n*（n-1）/2）;
  D = （double*）calloc（ndsizeof（double））;
  l=0;
  for （i=0;i<（n-1）;i++） {
    for （j=（i+1）;j      D[l] = （y[i]-y[j]）/（x[i]-x[j]）;
      l++;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4502  2010-12-07 16:47  MK\Mann_Kendall_parameters.c

     文件       1460  2010-09-08 17:00  MK\MK.FOR

     目录          0  2010-12-27 11:29  MK

----------- ---------  ---------- -----  ----

                 5962                    3