/*--------------------------------------------------------------------------------------------------------------------
  軟件功能：根據調和分析的結果，生成tecplot可讀入的文件，供tecplot畫同潮圖。

  輸入：輸入文件中grid.dat中存有緯度、經度、水深數據，xyap.dat是matlab調和分析的結果，其中兩列分別為振幅（m）和相位（度）

  輸出：同潮圖.dat

  作者：江文濱

  日期：2008-9-19
---------------------------------------------------------------------------------------------------------------------*/
#include "stdio.h"
#include "params.h"

void main()
{
	FILE *fpr,*fpw;
	int i,j,k,tt,mark[IM*JM];
	static double x[IM*JM],y[IM*JM],amp[IM*JM],phase[IM*JM];
	char temp[10];
	//  讀入xyap.dat
	fpr=fopen("xyap.dat","r");
	for (i=0;i<4;i++)
		fscanf(fpr,"%s",&temp);
	k=1;
	for (i=0;i<IM*JM;i++)
	{
		fscanf(fpr,"%lf %lf %lf %lf\n",&x[i],&y[i],&amp[i],&phase[i]);
		if(amp[i]!=-1000)
		{
			mark[i]=k;
			k++;
		}
	}
	tt=0;
	for(j=0;j<JM-1;j++)
		for(i=1;i<IM;i++)
		{
			if((amp[j*IM+i]!=-1000)&&(amp[j*IM+i-1]!=-1000)&&(amp[(j+1)*IM+i]!=-1000)&&(amp[(j+1)*IM+i-1]!=-1000))
			{
				if(phase[j*IM+i-1]<0)
					phase[j*IM+i-1]=360+phase[j*IM+i-1];
				tt++;
			}
		}
	//  輸出同潮圖文件
	fpw=fopen("cotidal.dat","w");
	fprintf(fpw,"Title=\"cotidal\"\n");
	fprintf(fpw,"Variables=\"x\",\"y\",\"a\",\"p\"\n");
	fprintf(fpw,"Zone n=%d,e=%d,f=fepoint,et=quadrilateral\n",k-1,tt);
	for (i=0;i<IM*JM;i++)
	{
		if(amp[i]!=-1000)
			fprintf(fpw,"%lf  %lf  %lf  %lf\n",x[i],y[i],amp[i]*100,phase[i]);
	}
	fprintf(fpw,"\n");
	for(j=0;j<JM-1;j++)
		for(i=1;i<IM;i++)
		{
			if((amp[j*IM+i]!=-1000)&&(amp[j*IM+i-1]!=-1000)&&(amp[(j+1)*IM+i]!=-1000)&&(amp[(j+1)*IM+i-1]!=-1000))
				fprintf(fpw,"%6d %6d %6d %6d\n",mark[j*IM+i-1],mark[j*IM+i],mark[(j+1)*IM+i],mark[(j+1)*IM+i-1]);
		}
	fclose(fpw);
	fclose(fpr);
}