#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fitsio.h"

#define LINEPIX 2100
#define LINEPIY 200

void printerror( int status);
double wl(int px){
  //return(1.373 + 1.102e-4*px - 1.314e-9*px*px - 5.402e-14*px*px*px);
  return(1.373 + (1.102e-4*px - 1.314e-9*px*px - 5.402e-14*px*px*px)*4.15);
}

double bb(int px){
  double w,dw,temp=1095+273,d=1240,s1=0.01,s2=1.5,tin=20,conv=2.1;

  w=wl(px);
  dw=wl(px+1)-w;
  return(1.191e11/w/w/w/w/w/(exp(14387.7/w/temp)-1)*dw*tin*(M_PI*s1/2*s1/2)*(M_PI*s2/2*s2/2)/d/d/(6.6260755e-27*2.99792458e+10/(w*1e-4))/conv);
}

int main(int argc, char *argv[])
{
  fitsfile *fp;
  int i,j,status,nfound,anynull;
  long fpixel;
  long naxes[2]={LINEPIY, LINEPIX},nbuffer;
  float nullval;
  static float buffer[LINEPIX],z[LINEPIY][LINEPIX];

  int px;
  float w,ws;

  if(argc!=3){
    fprintf(stderr,"Usage : %s input_file w_from(0.8815,1.1130,1.3815,1.5930)\n",argv[0]);
    exit(0);
  }
  sscanf(argv[2],"%f",&ws);

  fprintf(stderr,"Read %s  ... ",argv[1]);
  status=0;
  if(fits_open_file(&fp,argv[1], READONLY, &status))
    printerror(status);
  if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
    printerror(status);
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  if((naxes[0]>LINEPIX)||(naxes[1]>LINEPIY)){
    fprintf(stderr,"\nIllegal image format\n");
    exit(0);
  }
  fpixel = 1;
  nullval = 0;
  nbuffer = naxes[0];
  for(i=0;i<naxes[1];i++){
    if (fits_read_img(fp, TFLOAT, fpixel, nbuffer, &nullval,
                      buffer, &anynull, &status))
      printerror(status);
    for(j=0;j<naxes[0];j++) z[i][j]=buffer[j];
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  px=-5000;w=0;
  while(w<ws){
    w=wl(px++);
  }
  fprintf(stderr,"PIX=%d\n",px);
  for(j=0;j<naxes[0];j++){
    printf("%lf",wl(px+j));
    for(i=0;i<naxes[1];i++){
      printf(" %g",z[i][j]);
    }
    printf(" %lf",bb(px+j));
    printf("\n");
  }
  return(0);
}

void printerror( int status)
{
    /*****************************************************/
    /* Print out cfitsio error messages and exit program */
    /*****************************************************/

    char status_str[FLEN_STATUS], errmsg[FLEN_ERRMSG];

    if (status)
      fprintf(stderr, "\n*** Error occurred during program execution ***\n");

    fits_get_errstatus(status, status_str);   /* get the error description */
    fprintf(stderr, "\nstatus = %d: %s\n", status, status_str);

    /* get first message; null if stack is empty */
    if ( fits_read_errmsg(errmsg) )
    {
         fprintf(stderr, "\nError message stack:\n");
         fprintf(stderr, " %s\n", errmsg);

         while ( fits_read_errmsg(errmsg) )  /* get remaining messages */
             fprintf(stderr, " %s\n", errmsg);
    }

    exit( status );       /* terminate the program, returning error status */
}