/************************************************************************/ 
/* program to calculate the concentration from flux   */ 

/*   c = M * s
    M - response function  m x n
	s - estimate of s.  n x 1

	m: # of stations
	n: # of regions
*/ 

/* Ref: Deng, F. et al, 2007. Global monthly CO2 flux inversion with focus over North America. Tellus, 59B: 179-190 */ 

/* Written by G. Mo   June, 2012  */ 
/************************************************************************/ 

#include <stdio.h> 
#include <stdlib.h> 
//#include <math.h> 
#include <malloc.h> 
#include <string.h> 

int main()
{

/*  set matrix size here */
	int m=3066; 
	int n=3001; 
/*  ------------- */

	double **M,*s,*c;
	float v0,v1,v3; 
	
	int i,j,k;
	char in1_fn[80],in2_fn[80],outp_fn[180];
    	FILE *in1_ptr,*in2_ptr,*outp_ptr;

	M = (double **) malloc ((m+1)* sizeof(double));
		for (i=0; i<=m; i++)
		M[i] = (double *) malloc ((n+1)* sizeof(double));
	s = (double *) malloc (n* sizeof (double));
	c = (double *) malloc (m* sizeof (double));

/***   ***/
printf("reading M matrix ...\n");

    sprintf(in1_fn,"../response_73stn.bin");
    if ((in1_ptr=fopen(in1_fn, "rb")) == NULL)
    {
      printf("\n Unable to open file1,  exitting program ...\n\n");
	  exit(0);
    }

	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
		{
			fread(&v0,sizeof(float),1,in1_ptr);
			M[i][j]=v0;
		}
/***   ***/
printf("reading s matrix ...\n");

sprintf(in2_fn,"../flux_xd/firec13_flux_xd60m50r.txt");
//sprintf(in2_fn,"../flux_xd/ocnc13_flux_xd60m50r_pt2.txt");
//sprintf(in2_fn,"../flux_xd/bioc13_flux_xd60m50r_pt2.txt");
    if ((in2_ptr=fopen(in2_fn, "r")) == NULL)
    {
      printf("\n Unable to open file2,  exitting program ...\n\n");
	  exit(0);
    }
	for(i=0;i<n;i++)
	{
		fscanf(in2_ptr,"%f\n", &v1);
		s[i]=v1;
	}
	printf("s=%f \n",s[n-1]);;

/***    ***/
printf("Calculate M * s ...\n");

	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
		{
		  c[i]=0;
			for(k=0;k<n;k++)
				c[i] = c[i] + M[i][k]*s[k];
		}

/*  output inversed matrix*/
printf("output ...\n");

sprintf(outp_fn,"../d_xconcn4iterms/fxd_conctrn42x73f2c_fire.txt");
//sprintf(outp_fn,"../d_xconcn4iterms/fxd_conctrn42x73f2c_ocn_pt2.txt");
//sprintf(outp_fn,"../d_xconcn4iterms/fxd_conctrn42x73f2c_bio_pt2.txt");   
    if ((outp_ptr=fopen(outp_fn, "w")) == NULL)
    {
      printf("\n Unable to open output file,  exitting program ...\n\n");
	  exit(0);
    }

	for(i=0;i<m;i++)
		{
			v3=(float)c[i];
			fprintf(outp_ptr,"%f\n",v3);
		}
fclose(outp_ptr);
printf("done.\n");

return 1;
}   // end of main