#include <sys/types.h>
#include <unistd.h>

/* ============================================================================ */
/* module brightness_temp.c - convert observed IR radiance to temperature       */
/*                                                                              */
/* Written By: B. Franz, NASA/SIMBIOS, August 2003.                             */
/*                                                                              */
/* ============================================================================ */

#include "l12_proto.h"
#include "hdf_utils.h"

#define NTAB 100
#define NDET 10

static float  radtab  [NBANDSIR][NDET][NTAB];
static float  temptab [NBANDSIR][NDET][NTAB];
static float  radinc  [NBANDSIR][NDET];


/* ----------------------------------------------------------------------------------- */
/* read_bt_table() - reads the radiance to temperature HDF file                        */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
int read_bt_table(int32 sensorID)
{
    int32 sd_id;
    char *filedir;
    char filename[FILENAME_MAX];
    char sensor[20];

    /* Locate datafile using environment variable                   */
    if ((filedir = getenv("MSL12_DATA")) == NULL) {
        printf("-E- %s: MSL12_DATA env variable undefined.\n", __FILE__);
        return(0);
    }

    strcpy(sensor,&sensorName[sensorID][0]);
    lowercase(sensor);
    strcpy(filename, filedir);
    strcat(filename, "/");
    strcat(filename,sensor);
    strcat(filename,"/cal/");
    strcat(filename,"bt_");
    strcat(filename,sensor);
    strcat(filename,".hdf");

    /* Open the file */
    sd_id = SDstart(filename, DFACC_RDONLY);
    if(sd_id == FAIL){
        fprintf(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
        __FILE__,__LINE__,filename,DFACC_RDONLY);
        return(1);
    }

    READ_SDS_ID(sd_id,"Radiances"   ,radtab , 0,0,0,0, NBANDSIR,NDET,NTAB,1);
    READ_SDS_ID(sd_id,"Temperatures",temptab, 0,0,0,0, NBANDSIR,NDET,NTAB,1);

    printf("Radiance to brightness temperature loaded from %s\n",filename);

    SDend(sd_id);

    return(0);
}

/* ----------------------------------------------------------------------------------- */
/* radiance2bt() - interpolates brightness temps from table                .           */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
void radiance2bt(l1str *l1rec)
{
    static int firstCall = 1;

    long  id, ip, ib, ipb;
    long  i;
    float a;
    float rad;

    /* On first call, load the table */
    if (firstCall) {
        firstCall = 0;
        if (read_bt_table(l1rec->sensorID) != 0) {
 	    printf("Error loading brightness temperature table.\n");
	    exit(1);
	}
        for (ib=0; ib<NBANDSIR; ib++)
	    for (id=0; id<NDET; id++)
                radinc[ib][id] = radtab[ib][id][1] - radtab[ib][id][0];
    }

    id = l1rec->detnum;

    for (ip=0; ip<l1rec->npix; ip++)
        for (ib=0; ib<NBANDSIR; ib++) {

            ipb = ip*NBANDSIR+ib;
            rad = l1rec->Ltir[ipb]*10.0; /* radiance in W/m^2/nm/sr */

            if (rad <= radtab[ib][id][0])
		l1rec->Bt[ipb] = BT_LO;
            else if (rad >= radtab[ib][id][NTAB-1])
		l1rec->Bt[ipb] =  BT_HI;
            else {
                i = (long) ((rad - radtab[ib][id][0])/radinc[ib][id]);
                a = (rad-radtab[ib][id][i]) / (radtab[ib][id][i+1]-radtab[ib][id][i]);
                l1rec->Bt[ipb] = temptab[ib][id][i] + a*(temptab[ib][id][i+1]-temptab[ib][id][i]);
                l1rec->Bt[ipb] -= 273.15;
	    }
	}
}


/* ----------------------------------------------------------------------------------- */
/* get_sst_modist() - compute long wave sea surface temperature for MODIS/Terra        */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
float get_sst_modist(float Bt11, float Bt12, float senz, float sstref)
{
    static float pi = 3.141592654;
    static float dBtlo = 0.5;
    static float dBthi = 0.9;
    static float a[2][4] = {{1.0520,0.9840,0.1300,1.860},
                            {1.8860,0.9380,0.1280,1.094}};
    float dBt = Bt11 - Bt12;
    float mu  = cos(senz*pi/180.0);
    float sstlo, ssthi, sst;

    if (dBt <= dBtlo)
        sst = a[0][0] + a[0][1]*Bt11 + a[0][2]*dBt*sstref + a[0][3]*dBt*(1.0/mu-1.0);
    else if (dBt >= dBthi)
        sst = a[1][0] + a[1][1]*Bt11 + a[1][2]*dBt*sstref + a[1][3]*dBt*(1.0/mu-1.0);
    else {
        sstlo = a[0][0] + a[0][1]*Bt11 + a[0][2]*dBt*sstref + a[0][3]*dBt*(1.0/mu-1.0);
        ssthi = a[1][0] + a[1][1]*Bt11 + a[1][2]*dBt*sstref + a[1][3]*dBt*(1.0/mu-1.0);
        sst = sstlo + (dBt-dBtlo)/(dBthi-dBtlo)*(ssthi-sstlo);
    }
 
    return(sst);
}

/* ----------------------------------------------------------------------------------- */
/* get_sst_modisa() - compute long wave sea surface temperature for MODIS/Aqua         */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
float get_sst_modisa(float Bt11, float Bt12, float senz, float sstref)
{
    static float pi = 3.141592654;
    static float dBtlo = 0.5;
    static float dBthi = 0.9;
    static float a[2][4] = {{1.1520,0.9600,0.1510,2.0210},
                            {2.1330,0.9260,0.1250,1.1980}};
    float dBt = Bt11 - Bt12;
    float mu  = cos(senz*pi/180.0);
    float sstlo, ssthi, sst;

    if (dBt <= dBtlo)
        sst = a[0][0] + a[0][1]*Bt11 + a[0][2]*dBt*sstref + a[0][3]*dBt*(1.0/mu-1.0);
    else if (dBt >= dBthi)
        sst = a[1][0] + a[1][1]*Bt11 + a[1][2]*dBt*sstref + a[1][3]*dBt*(1.0/mu-1.0);
    else {
        sstlo = a[0][0] + a[0][1]*Bt11 + a[0][2]*dBt*sstref + a[0][3]*dBt*(1.0/mu-1.0);
        ssthi = a[1][0] + a[1][1]*Bt11 + a[1][2]*dBt*sstref + a[1][3]*dBt*(1.0/mu-1.0);
        sst = sstlo + (dBt-dBtlo)/(dBthi-dBtlo)*(ssthi-sstlo);
    }
 
    return(sst);
}


/* ----------------------------------------------------------------------------------- */
/* get_sst() - compute sea surface temperature                             .           */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
float comp_sst(l2str *l2rec, long ip)
{
    static float Btmin   = -3.0;
    static float Btmax   = 35.0;
    static float SSTdiff =  6.0;
    static float dBtmin  = -1.1;
    static float dBtmax  = 10.0;  

    float sst    = SST_BAD;
    float sstref = l2rec->sstref[ip];
    long  ipb    = ip*NBANDSIR;
    float Bt11   = l2rec->Bt[ipb];
    float Bt12   = l2rec->Bt[ipb+1];
    float senz   = l2rec->senz[ip];
    float dBt    = Bt11-Bt12;

    /* if brights are flagged as bogus, flag SST and return */
    if (Bt11 < BT_LO+0.1 || Bt11 > BT_HI-0.1 || 
        Bt12 < BT_LO+0.1 || Bt12 > BT_HI-0.1) {
        l2rec->flags[ip] |= SSTFAIL;
        return (sst);
    }

    switch (l2rec->sensorID) {
        case MODIST:
            sst = get_sst_modist(Bt11,Bt12,senz,sstref);
            break;
        case MODISA:
            sst = get_sst_modisa(Bt11,Bt12,senz,sstref);
            break;
        default:
            l2rec->flags[ip] |= SSTFAIL;
            return (sst);
    }

    if ( Bt11 <  Btmin ||  Bt11 >  Btmax || 
         Bt12 <  Btmin ||  Bt12 >  Btmax ||
         dBt  < dBtmin ||  dBt  > dBtmax ||
         abs(sst-sstref) > SSTdiff ) {

        l2rec->flags[ip] |= SSTWARN;
    }

    return (sst); 
}



/* ----------------------------------------------------------------------------------- */
/* get_oisst() - read and interpolate flat binary Reynolds OIV2 SST                    */
/*                                                                                     */
/* The files were written in IEEE binary (big-endian). Each file contains four FORTRAN */
/* records described as follows:                                                       */
/*                                                                                     */
/*    rec 1: date and version number        (8 4-byte integer words)                   */
/*    rec 2: gridded sst values in degC     (360*180 4-byte real words)                */
/*    rec 3: normalized error variance      (360*180 4-byte real words)                */
/*    rec 4: gridded ice concentration      (360*180 1-byte integer words)             */
/*                                                                                     */
/* B. Franz, SAIC, May 2004.                                                           */
/* ----------------------------------------------------------------------------------- */

#define OINX 360
#define OINY 180

float get_oisst(char *sstfile, float lon, float lat)
{
    static int   firstCall = 1;
    static int   nx = OINX;
    static int   ny = OINY;
    static float dx = 360.0/OINX;
    static float dy = 180.0/OINY;
    static float sstref[OINY+2][OINX+2];

    float sst  = SST_BAD;
    int   i,j,ii;
    float xx,yy;
    float t,u;

    if (firstCall) {

        FILE *fp = NULL;
        float ssttmp[OINY][OINX];
        long syear,smon,sday;
        long eyear,emon,eday;
        long ndays,version;

        firstCall = 0;

        if ((fp = fopen(sstfile,"r")) == NULL) {
            printf("Error opening SST reference file %s for reading.\n",sstfile);
            exit(1);
        }

        if (fseek(fp,4,SEEK_SET) < 0) {
            printf("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}          
        if (fread(&syear   ,sizeof(long),1,fp) != 1) {
            printf("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}
        fread(&smon    ,sizeof(long),1,fp);
        fread(&sday    ,sizeof(long),1,fp);
        fread(&eyear   ,sizeof(long),1,fp);
        fread(&emon    ,sizeof(long),1,fp);
        fread(&eday    ,sizeof(long),1,fp);
        fread(&ndays   ,sizeof(long),1,fp);
        fread(&version ,sizeof(long),1,fp);
        fseek(fp,4,SEEK_CUR);

#ifdef linux
        swapc_bytes(&syear    ,sizeof(long),1);
        swapc_bytes(&smon     ,sizeof(long),1);
        swapc_bytes(&sday     ,sizeof(long),1);
        swapc_bytes(&eyear    ,sizeof(long),1);
        swapc_bytes(&emon     ,sizeof(long),1);
        swapc_bytes(&eday     ,sizeof(long),1);
        swapc_bytes(&ndays    ,sizeof(long),1);
        swapc_bytes(&version  ,sizeof(long),1);
#endif

        printf("Loading OI Reynolds SST reference from %s\n",sstfile);
        printf("  file start date (y/m/d): %d / %d / %d\n",syear,smon,sday);
        printf("  file end   date (y/m/d): %d / %d / %d\n",eyear,emon,eday);
        printf("  days composited: %d\n",ndays);
        printf("  file version: %d\n",version);
        printf("\n");

        if (fseek(fp,4,SEEK_CUR) < 0) {
            printf("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}          
        if (fread(&ssttmp[0][0],sizeof(float),nx*ny,fp) != nx*ny) {
            printf("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}
        fclose(fp);

#ifdef linux
        swapc_bytes(&ssttmp[0][0],4,nx*ny);
#endif

        /* rotate 180-deg and add wrapping border to simplify interpolation */
        /* new grid is -180.5,180.5 in i=0,361 and -90.5,90.5 in j=0,181    */

        for (j=0; j<ny; j++) {
 	    for (i=0; i<nx; i++) {
	        ii = (i < nx/2) ?  i+nx/2 : i-nx/2;
                sstref[j+1][ii+1] = ssttmp[j][i];
 	    }
            sstref[j+1][0]    = sstref[j+1][nx];
            sstref[j+1][nx+1] = sstref[j+1][1];
        }
        for (i=0; i<nx+2; i++) {
            sstref[0]   [i] = sstref[1][i];
            sstref[ny+1][i] = sstref[ny][i];
        }

	/*
        fp = fopen("/tmp/sst.dat","w");
        fwrite(&sstref[0][0],4,(nx+2)*(ny+2),fp);
        fclose(fp);
        printf("writing test sst and exiting\n");
        exit(1);
	*/
    }


    /* locate LL position within reference grid */
    i = MAX(MIN((int) (lon+180.0+dx/2)/dx,OINX+1),0);
    j = MAX(MIN((int) (lat+ 90.0+dy/2)/dy,OINY+1),0);

    /* compute longitude and latitude of that grid element */
    xx = i*dx - 180.0 - dx/2;
    yy = j*dy -  90.0 - dy/2;

    /* best match to miami */
    /*
    i = MAX(MIN((int) (lon+180.0+dx)/dx,OINX+1),0);
    j = MAX(MIN((int) (lat+ 90.0   )/dy,OINY+1),0);
    xx = i*dx - 180.0 - dx;
    yy = j*dy -  90.0;
    */

    /* bilinearly interpolate */
    t = (lon - xx)/dx;
    u = (lat - yy)/dy;

    sst = (1-t)*(1-u) * sstref[j  ][i  ]
        + t*(1-u)     * sstref[j  ][i+1]
        + t*u         * sstref[j+1][i+1]
        + (1-t)*u     * sstref[j+1][i  ];

    /*
    sst = sstref[j][i];
    */

    return(sst);    
}



/* ----------------------------------------------------------------------------------- */
/* get_sstref() - retrieves reference sea surface temperature              .           */
/*                                                                                     */
/* B. Franz, SAIC, May 2004.                                                           */
/* ----------------------------------------------------------------------------------- */
#define PATHCLIM 1
#define OISSTBIN 2
#include "smi_climatology.h"
float get_sstref(char *sstfile, float lon, float lat, int day)
{
    static long format = -1;

    float sst    = SST_BAD;
    long  sd_id;

    if (format < 0) {

        /* Does the file exist? */
        if (access(sstfile, F_OK) || access(sstfile, R_OK)) {
            printf("-E- %s: SST input file '%s' does not exist or cannot open.\n",
                   __FILE__, sstfile);
            exit(1);
        }

        /* What is it? */
        sd_id = SDstart(sstfile, DFACC_RDONLY);
        if (sd_id != FAIL) {

            /* Format is HDF,  Assuming Pathfindaer climatology. */
  	    format = PATHCLIM;
            if (smi_climatology_init(sstfile,day,SST) != 0) {
                printf("-E- %s : Unable to initialize SST file\n",__FILE__);
                printf("%s %d\n",sstfile,day);
                exit(1);
	    }             

        } else {

	    /* Format is not HDF.  Assuming OISST Binary. */
  	    format = OISSTBIN;
        }

    }

    switch (format) {
      case PATHCLIM:
        sst = smi_climatology(lon,lat,SST);
        break;
      case OISSTBIN:
        sst = get_oisst(sstfile,lon,lat);
        break;
      default:
        printf("-E- %s: unknown SST input file format for %s.\n",__FILE__, sstfile);
        break;
    }
    return(sst);
}