/*--------------------------------------------------------------------------------------------
  File name: SMOKE2AREA.c

  Purpose:
  	Read Goes East McIDAS AREA file as template and creat one for Smoke AOD
  
  Adapted from AOD2AREA.c which was written By Zedong Zhang on 2/18/2004
  Edited by Laura Kowal to work for SMOKE data April 2008
  Modified: By J.Yang in May 2008

  Usage:
    	SMOKE2AREA <Smoke binary file name>

  Notes:
       The Smoke binary data values are as following:
     	Smoke AOD:                   > 0.  (Maximum is 1.517, we use 1.6 for scaling)
     	Cloud:                       -1.
     	No fire:                     -2.
     	Nonsmoke aerosols:           -3.
     	No data:                     -9999.

       In McIDAS AREA file the brightnesss values are as following:
	no data:  	   	0
	cloud:    	   	1
	no fire:  	   	2
	nonsmoke aerosol:  	3 
        Smoke AOD:           10-255 


  Last Modified: 05/20/2008
--------------------------------------------------------------------------------------------*/

#include <stdio.h>
#include <string.h>
#include <time.h>
#define WORD unsigned int
#define BYTE unsigned char
#define SHORT unsigned short

#define MX 2126
#define MY 880
#define DX 2
#define DY 14
#define Image_WEST_X 2130

	
struct Directoryblock {
  WORD relp;     //1. relative position of the image object in the ADDE dataset
  WORD imgtype;  //2. image type=4
  WORD sennum;   //3. SSEC sensor source number
  WORD yyyddd;   //4. nominal year and Julian day of the image, yyyddd
  WORD hhmm;     //5. nominal time of the image, hhmmss
  WORD lineC;    //6. upper-left image line coordinate
  WORD eleC;     //7. upper-left image element coordinate
  WORD res1;     //8. reserved;:w
  WORD numline;  //9. number of line in the image
  WORD numdata;  //10. number of data points per line
  WORD numbyte;  //11. number of bytes per data point
  WORD lineresol;//12. line resolution
  WORD eleres;   //13. element resolution
  WORD numband;  //14. number of spectral bands;
  WORD lenprefix;//15. length of the line prefix
  WORD projnum;  //16. SSEC file creation project number
  WORD filedate; //17. file creation year and Julian day, yyyddd
  WORD filetime; //18. file creation time, hhmmss
  WORD bandmap;  //19. spectral band map
  WORD ID;       //20. image ID number
  char res2[16]; //21-24. resrved for sensor-specific data
  char memo[32]; //25-32. memo field; 32 ASCII characters
  WORD res3;     //33. reserved
  WORD datablock;//34. byte offset to the start of the data block
  WORD navblock; //35. byte offset to the start of the navigation block
  WORD valcode;  //36. validity code;
  char PDL[32];  //37-44. Program Date Load; used for pre-GOES-8 satellites
  WORD band8;    //45. source of band 8; used for GOES AA processing
  WORD imgdate;  //46. actual image start year and Julian day, yyyddd
  WORD imgtime;  //47. actual image start time, hhmmss; in milliseconds for POES data
  WORD imgscan;  //48. actual image start scan
  WORD lendoc;   //49. length of the prefix documentation
  WORD lencali;  //50. length of the prefix calibration
  WORD lenband;  //51. length of the prefix band list
  char stype[4]; //52. source type; satellite-specific (ASCII)
  char ctype[4]; //53. calibration type; satellite-specific (ASCII)
  char res4[24]; //54-59. reserved
  WORD supblock; //60. byte offset to the supplemental block
  WORD numsup;   //61. number of bytes in the supplemental block
  char res5[4];  //62. reserved
  WORD caliblock;//63. byte offset to the start of the calibration block
  WORD numcom;   //64. number of comment cards
}dheader;

void usage(void) {
	printf("Usage: SMOKE2AREA SMOKE_file\n");
}

int searchSlash(char *s)
{
   int i,p;

   p=0;
   for(i=0;i<strlen(s);i++)
      if(s[i]=='/') p=i+1;

   return p;
}

// Julian day conversion function

int getJulianDay(int, int, int);

int getJulianDay(int year, int month, int day) {

   int jday;
   int noleapyear[12] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
   int leapyear[12]   = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335};

   if (((year % 4 == 0) && (year % 100 != 0)) || (year % 400 == 0))
       jday = leapyear[month-1] + day;
   else
       jday = noleapyear[month-1] + day;

   return jday;
}
// main program

int main(int argc,char *argv[])
{
   FILE *fp,*ft;
   WORD i, j, k;  
   BYTE imagedata,zero=0;
   WORD Offset;   
   float data[MX][MY];
   char accmask;
   double x;
   time_t timer;  
   struct tm mytime; 
   int img_x;
   int left_x;
   int img_y;


   // minimum count for the '0' Smoke concentration
   int minCnt = 10;   

   if(argc<2) {
	usage();
   	exit(1);
   }

   printf("Read AREA file\n");
   //printf("***************************\n");
/*
   fp=fopen("/data/gasp/work/AREA_template","rb");

   if(fp==NULL) {
    printf("Error: AREA_template doesn't exist!\n");
    exit(1);
  }
  
  if(fread(&dheader,sizeof(dheader),1,fp)!=1) {
    printf("Error to read this file!\n");
    fclose(fp);
    exit(1);
  }

  fclose(fp);

  system("cp /data/gasp/work/AREA_template AOD_AREA");
*/
  fp=fopen("/data/SMOKE/work/SMOKE_AREA","rb+");
   
  if(fread(&dheader,sizeof(dheader),1,fp)!=1) {
    printf("Error to read this file!\n");
    fclose(fp);
    exit(1);
  }

//  printf("Bands:%d; Number Byte per Data:%d\n", dheader.numband,dheader.numbyte);
//  printf("Offset of Data Block:%d; Number of Lines:%d; Number of elements:%d;\n",dheader.datablock,dheader.numline,dheader.numdata);


  //Read date from this AREA file
/*  if(fseek(fp, dheader.datablock, SEEK_SET)!=0) {
    printf("Error to read this file!\n");
    fclose(fp);
    exit(1);
  }
*/

   int yyyymmdd;
   int julian;
   int yyyy;
   int mm;
   int dd;
   int yyy;
   int hhmm;
   int imgdate;
   int imgtime;
   


  julian=0;
//  * = &cccddyyy;

  yyyymmdd=0;
  j=1;
 // printf("%d\n",searchSlash(argv[1]));
  for(i=0;i<8;i++) {
     yyyymmdd+=(argv[1][searchSlash(argv[1])+11-i]-'0')*j;
      //printf("%d\n", yyyymmdd);

  j=j*10;
  }
  // year
  j=1;
  yyyy=0;
  for(i=0;i<4;i++) {
     yyyy+=(argv[1][searchSlash(argv[1])+7-i]-'0')*j;
      //printf("%d\n", yyyy);

  j=j*10;
  }
 
//  year = &yyyy;

//printf("%d\n", year);

  //month

  j=1;
  mm=0;
     for(i=0;i<2;i++) {
     mm+=(argv[1][searchSlash(argv[1])+9-i]-'0')*j;
      //printf("%d\n", mm);

  j=j*10;
  }

//  month = &mm;

  //day
    j=1;
    dd=0;
     for(i=0;i<2;i++) {
     dd+=(argv[1][searchSlash(argv[1])+11-i]-'0')*j;
      //printf("%d\n", dd);

   j=j*10;
   }

//   day = &dd;
   //year for julian date

       j=1;
       yyy=0;
    for(i=0;i<3;i++) {
      yyy+=(argv[1][searchSlash(argv[1])+7-i]-'0')*j;
      //printf("%d\n", yyy);

    j=j*10;
    }
      yyy+=100;
   //printf("%d\n", yyy);
  
  // hhmm

  
  hhmm=0;
  j=1;
   for(i=0;i<4;i++) {
      hhmm+=(argv[1][searchSlash(argv[1])+15-i]-'0')*j;
      j=j*10;
   }

  hhmm *= 100; 
  //printf("%d\n",hhmm);
 

// Calculate Julian Day
 
  int year=yyyy, month=mm, day=dd;  
  int JulianDay;
    

  JulianDay = getJulianDay(year, month, day);
  printf("Julian day %d for Year %d Month %d and Day %d\n",JulianDay,year, month, day);

// display time

  int yyyddd;

    yyyddd = (yyy*1000)+JulianDay;

printf("McIDAS Date:%d Time:%d\n", yyyddd,hhmm);
  imgdate=yyyddd;
  dheader.imgdate=yyyddd;
  dheader.yyyddd=yyyddd;
  dheader.hhmm=hhmm;
  dheader.imgtime=hhmm;
  imgtime=dheader.hhmm;

//printf("%d\n", dheader.hhmm);
//printf("%d\n", imgtime);

  timer=time(NULL);
  mytime=*gmtime(&timer);


  printf("Current Date: %d %d\n",(mytime.tm_year+1900)*1000+mytime.tm_yday+1,mytime.tm_hour*10000+mytime.tm_min*100+mytime.tm_sec);


//  dheader.bandmap=0;     

  dheader.filedate=mytime.tm_year*1000+mytime.tm_yday+1;
  dheader.filetime=mytime.tm_hour*10000+mytime.tm_min*100+mytime.tm_sec;


  strcpy(dheader.memo,"SMOKE AOD  ");  
  //memcpy (dheader.stype, "PRD ", 4);  
  //memcpy (dheader.ctype, "BRIT", 4);  

  fseek(fp,0,SEEK_SET);
  fwrite(&dheader,sizeof(dheader),1,fp);           


  ft=fopen(argv[1],"rb");
  for(j=0;j<MY;j++)
  for(i=0;i<MX;i++)
     //fread(&data[i][j],4,1,ft);
     fread(&data[i][j],sizeof(float),1,ft);

                                                                              
//  fseek(ft,4*MX*MY,SEEK_CUR);
                                                                               
//  for(j=0;j<MY;j++)
//  for(i=0;i<MX;i++) {
//     fread(&accmask,1,1,ft);
//     data[i][j]=data[i][j]*accmask;
//  }
  
  fclose(ft);



  if(fseek(fp, dheader.datablock, SEEK_SET)!=0) {
    printf("Error to read this file!\n");
    fclose(fp);
    exit(1);
  }
  img_y = MY + DY;
  left_x =Image_WEST_X - DX -MX;

  //Write to Image file
  // x=255/1.6;
  x = (255 - minCnt)/1.6;
//**  for(j=0;j<822;j++) {
  for(j=0;j<img_y;j++) {
    if(fseek(fp, dheader.lenprefix, SEEK_CUR)!=0) {
      printf("Error to read this file!\n");
      fclose(fp);
      exit(1);
    }
    if (j < DY) {
       for(i=0;i<Image_WEST_X;i++) {
          fwrite(&zero,sizeof(BYTE),1,fp);
       }
     }
    else {

    for(i=0;i<DX;i++)
         fwrite(&zero,sizeof(BYTE),1,fp);

   }

//    fseek(fp,7*sizeof(BYTE), SEEK_CUR);
//    for(i=0;i<7;i++)
//       fwrite(&zero,sizeof(BYTE),1,fp);

    for(i=0;i<MX;i++) {
//      x=data[i][822-j]/1000.;
//      if(data[i][822-j]>1300) x=0;
//       data[i][822-j] = data[i][822-j] / 1000000;
       //if((data[i][822-j]>0) && ((i % 10) == 0) )
       //{
       //    printf("data value = %f\n", data[i][822-j]);
       //}

      if(data[i][MY-j-DY]<0)
      {
         if(data[i][MY-j-DY] == -1.0)
         {
            imagedata = 1;
         }
         else if (data[i][MY-j-DY] == -2.0)
         {
            imagedata = 2;
           // if((i % 10) == 0)
           // {
           //   printf("-2 value = %f, count value is %d\n", data[i][822-j], imagedata);
           // }

         }
         else if (data[i][MY-j-DY] == -3.0)
         {
            imagedata = 3;
         }
         else
         {
            imagedata = 0;
            /* for debug
            if((i % 50) == 0)
            {
              printf("Other negative value = %f, count value is %d\n", data[i][822-j], imagedata);
            } */

         }
      }   
      else if (data[i][MY-j-DY]>1.6)
      {
 	 imagedata= 255;
      }
      else 
      {
        // plus 0.5 for round the value to the nearest integer 
        imagedata=(BYTE)(x*data[i][MY-j] + minCnt + 0.5);  

        /* for debug
        if((i % 50) == 0)
        {
           printf("Smoke data value = %f, count value is %d\n", data[i][822-j], imagedata);
        }
        */  

      } 

      fwrite(&imagedata,sizeof(BYTE),1,fp);
    }
//    fseek(fp,121*sizeof(BYTE), SEEK_CUR);      
//    for(i=0;i<121;i++)
//       fwrite(&zero,sizeof(BYTE),1,fp);
    

  }  

//  for(j=822;j<883;j++) 
//    for(i=0;i<2128;i++)
//       fwrite(&zero,sizeof(BYTE),1,fp);
   for(j=MY;j<dheader.numline;j++)
    for(i=0;i<dheader.numdata;i++)
       fwrite(&zero,sizeof(BYTE),1,fp);

  

  fclose(fp);

/*  fp=fopen("AOD_ASCII.dat","w");
  for(j=0;j<MY;j++) {
     for(i=0;i<MX;i++)
        fprintf(fp,"%d ",data[i][j]);
     fprintf(fp,"\n");
  }
  fclose(fp);
*/

//printf("%d\n", imgtime);

  exit(0);
}