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

#define LINEPIX 2100
#define LINEPIY 1800
#define EFFICIENCY 2.5

void printerror( int status );

int main(int argc, char *argv[])
{
  fitsfile *fp;
  int i,j,status,nfound,anynull,bitpix=-32;
  long fpixel, nelements,naxis=2;
  long naxes[2]={LINEPIX, LINEPIY},nbuffer;
  float nullval;
  static float buffer[LINEPIX],z[LINEPIY][LINEPIX],w[LINEPIY/9][LINEPIX];
  static float zo[LINEPIY][LINEPIX],wt[LINEPIY],tr[LINEPIY/9][LINEPIX],ab[LINEPIX],atm[LINEPIX];
  static float no[LINEPIY][LINEPIX],n[LINEPIY][LINEPIX];
  float wlmin,wlmax,dw,w1,w2,ws,we,ww,abt,abt2,atmt;
  float w0min[7]={1.67, 1.59, 1.51, 1.22, 1.11, 1.00, 1.28};
  float w0max[7]={1.73, 1.65, 1.57, 1.28, 1.17, 1.06, 1.34};
  float w1min[7]={1.64, 1.56, 1.48, 1.19, 1.08, 0.97, 1.16};
  float w1max[7]={1.66, 1.58, 1.50, 1.21, 1.10, 0.99, 1.26};
  float w2min[7]={1.74, 1.66, 1.58, 1.29, 1.18, 1.07, 1.50};
  float w2max[7]={1.76, 1.68, 1.60, 1.31, 1.20, 1.09, 1.60};
  float ff0[7]  ={2.00, 2.00, 2.00, 2.00, 2.30, 3.50, 1.00};
  float ff1[7]  ={0.91, 0.86, 1.10, 1.16, 1.85, 1.75, 1.26};
  float rf0[7]  ={0.92, 0.85, 0.89, 0.85, 0.80, 0.74, 0.83};
  float rf1[7]  ={1.00, 1.33, 1.25, 0.80, 0.96, 0.95, 1.00};
  float wei[9]  ={0,1,1,1,1,1,1,1,0},tmpval,bx[5]={0.1,0.2,0.4,0.2,0.1},af[30],t1,af3,af5,daf=0.;
  int nj,jw,iw,std,mode,ot,raw=0,t2,naxes0,naxes1,irs;
  double f1,f2,f0,a,af0,afm,expmin,trs1,trs2,fac,fac1,fac2,fac3,am=1.2,ff1t;
  FILE *fpd,*fps;
  char buf[256],fname[128],fname2[128],dmy[16];

  if(argc>5 && (strcmp(argv[5],"lr")==0 || strcmp(argv[5],"js")==0 || strcmp(argv[5],"jl")==0 || strcmp(argv[5],"hs")==0 || strcmp(argv[5],"hl")==0)) ot=1;
  else ot=0;

  if(argc!=5 && argc!=7+ot && argc!=8+ot){
    fprintf(stderr,"Usage : %s input_file WL_map cont_image output_file [exp_min #STAR [Airmass or Type or NONE]]\n",argv[0]);
    exit(0);
  }
  if(argc==5) {raw=1;std=100;}
  else{
    sscanf(argv[5+ot],"%lf",&expmin);
    if(strcmp(argv[6+ot],"RAW")==0 || strcmp(argv[6+ot],"raw")==0) {raw=1;std=100;}
    else sscanf(argv[6+ot],"%d",&std);
    if(argc==8+ot && (strcmp(argv[7+ot],"RAW")==0 || strcmp(argv[7+ot],"raw")==0)) raw=1;
  }
  if(std<1 || std>LINEPIY/9){
    fprintf(stderr,"#STAR must be between 1 and %d\n",LINEPIY/9);
    exit(0);
  }
  sscanf(argv[1],"%s",fname);
  fname[strlen(argv[1])-5]='\0';
  sprintf(fname2,"%s0.dat",fname);
  sprintf(fname,"%s0_nosq.fits",fname);

  if(strncmp(argv[1],"irs2",4)==0) irs=2;
  else irs=1;

  if((fpd=fopen("extract.dat","r"))==NULL){
    fprintf(stderr,"Can not open file: extract.dat\n Use default value.\n");
  }else{
    fprintf(stderr,"Open kernel file: extract.dat\n");
    for(i=0;i<9;i++){
      if(fscanf(fpd,"%f",&wei[i])!=1){
        fprintf(stderr,"Wrong format\n");
        exit(1);
      }
    }
   fclose(fpd);
  }
  if((fpd=fopen("conv_w_fnu.dat","r"))!=NULL){
    fprintf(stderr,"Open color correction file: conv_w_fnu.dat\n");
    for(i=0;i<7;i++){
      if(fscanf(fpd,"%f",&ff1[i])!=1){
        fprintf(stderr,"Wrong format\n");
        exit(1);
      }
    }
    for(i=0;i<7;i++){
      if(fscanf(fpd,"%f",&rf1[i])!=1){
        fprintf(stderr,"Wrong format\n");
        exit(1);
      }
    }
    fclose(fpd);
  }
  
  ww=0;
  fprintf(stderr,"Read %s  ... ",argv[3]);
  status=0;
  if(fits_open_file(&fp,argv[3], 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){
    fprintf(stderr,"\nNumber of columns must be < %d\n",LINEPIX);
    exit(0);
  }
  naxes0=naxes[0];
  if(naxes[1]!=LINEPIY/9){
    fprintf(stderr,"\nNumber of lines must be %d\n",LINEPIY/9);
    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);
    wt[i]=0;
    for(j=0;j<naxes[0];j++) {wt[i]+=buffer[j];ww+=buffer[j]/naxes[1];}
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  for(i=0;i<naxes[1];i++){
    wt[i]/=ww;
    //fprintf(stderr,"%d\t%g\n",i,wt[i]);getchar();
  }

  if((fpd=fopen(fname2,"r"))==NULL){
    fprintf(stderr,"Can not open file: %s , Set X-shift=0\n",fname2);
    ww=0.;
  }else{
    fprintf(stderr,"Read %s ... ",fname2);
    fscanf(fpd,"%f\n",&ww);
    fprintf(stderr,"X-shift=%g\n",ww);
    fclose(fpd);
  }
  jw=(int)floor(ww);
  wlmin=1.8;wlmax=0.9;
  fprintf(stderr,"Read %s  ... ",argv[2]);
  status=0;
  if(fits_open_file(&fp,argv[2], 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]!=naxes0){
    fprintf(stderr,"\nNumber of columns is different from that of %s (%d)\n",argv[3],naxes0);
    exit(0);
  }
  if(naxes[1]!=LINEPIY/9){
    fprintf(stderr,"\nNumber of lines must be %d\n",LINEPIY/9);
    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++){
      if(j+jw<0) w[i][j]=buffer[1]*(ww+j)+buffer[0]*(1-ww-j);
      else if(j+jw+1>=naxes[0]) w[i][j]=buffer[naxes[0]-1]*(ww-naxes[0]+2+j)+buffer[naxes[0]-2]*(1-ww+naxes[0]-2-j);
      else w[i][j]=buffer[j+jw+1]*(ww-jw)+buffer[j+jw]*(1-ww+jw);
      if(wlmin>w[i][j]) wlmin=w[i][j];
      if(wlmax<w[i][j]) wlmax=w[i][j];
      //fprintf(stderr,"%d %d %g %g %g %g\n",i,j,buffer[j],w[i][j],wlmin,wlmax);getchar();
    }
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");
  naxes0=naxes[0];naxes1=naxes[1];
  fprintf(stderr,"wlmin=%g wlmax=%g\n",wlmin,wlmax);

  mode=0;
  for(i=1;i<8;i++){
    if(wlmin<w1min[i-1] && w2max[i-1]<wlmax){
      if(mode!=0){
	if(fabs(wlmax+wlmin-w2max[i-1]-w1min[i-1])<fabs(wlmax+wlmin-w2max[mode-1]-w1min[mode-1])) mode=i;
      }else mode=i;
    }
  }
  if(mode==0){
    fprintf(stderr,"Observed wavelength (%g - %g) must contain one of the following range.\n",wlmin,wlmax);
    fprintf(stderr,"     J-short : %g - %g\n",w1min[5],w2max[5]);
    fprintf(stderr,"     J-middle: %g - %g\n",w1min[4],w2max[4]);
    fprintf(stderr,"     J-long  : %g - %g\n",w1min[3],w2max[3]);
    fprintf(stderr,"     H-short : %g - %g\n",w1min[2],w2max[2]);
    fprintf(stderr,"     H-middle: %g - %g\n",w1min[1],w2max[1]);
    fprintf(stderr,"     H-long  : %g - %g\n",w1min[0],w2max[0]);
    fprintf(stderr,"     Low res : %g - %g\n",w1min[6],w2max[6]);
    exit(1);
  }
  if(mode==7){wlmin=0.9;wlmax=1.8;dw=0.0005;}
  else{
    wlmin=floor((wlmin+wlmax)/2*1000+0.5)/1000-0.12;
    wlmax=wlmin+0.24;
    dw=0.000125;
    if(wlmin<0.9){wlmin=0.9;wlmax=1.14;}
    if(wlmax>1.8){wlmin=1.56;wlmax=1.80;}
  }
  if(mode==1) fprintf(stderr,"H-long :");
  else if(mode==2) fprintf(stderr,"H-middle :");
  else if(mode==3) fprintf(stderr,"H-short :");
  else if(mode==4) fprintf(stderr,"J-long :");
  else if(mode==5) fprintf(stderr,"J-middle :");
  else if(mode==6) fprintf(stderr,"J-short :");
  else fprintf(stderr,"Low res :");
  fprintf(stderr," wlmin=%g wlmax=%g dw=%g\n",wlmin,wlmax,dw);

  if(raw!=1){
    fprintf(stderr,"Read %s  ... ",fname);
    status=0;
    if(fits_open_file(&fp, fname, READONLY, &status)){
      printerror(status);
      //fprintf(stderr,"failed!\n");
      //fprintf(stderr,"Input noise^2 level is assumed to be 1.\n");
      //for(i=0;i<LINEPIY;i++){
      //for(j=0;j<naxes[0];j++) n[i][j]=1.;
      //}
    }else{
      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]!=naxes0){
	fprintf(stderr,"\nNumber of columns is different from that of %s(%d)\n",argv[3],naxes0);
	fprintf(stderr,"Input noise^2 level is assumed to be 1.\n");
	for(i=0;i<LINEPIY;i++){
	  for(j=0;j<naxes[0];j++) n[i][j]=1.;
	}      
      }else if(naxes[1]!=LINEPIY&&naxes[1]!=LINEPIY/9){
	fprintf(stderr,"\nNumber of lines must be %d or %d\n",LINEPIY,LINEPIY/9);
	fprintf(stderr,"Input noise^2 level is assumed to be 1.\n");
	for(i=0;i<LINEPIY;i++){
	  for(j=0;j<naxes[0];j++) n[i][j]=1.;
	}      
      }else{
	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++) n[i][j]=buffer[j];
	  fpixel+=naxes[0];
	}
      }
      if(fits_close_file(fp, &status))
	printerror(status);
      fprintf(stderr,"OK\n");
    }
  }

  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]!=naxes0){
    fprintf(stderr,"\nNumber of columns is different from that of %s(%d)\n",argv[3],naxes0);
    exit(0);
  }
  if(naxes[1]!=LINEPIY&&naxes[1]!=LINEPIY/9){
    fprintf(stderr,"\nNumber of lines must be %d or %d\n",LINEPIY,LINEPIY/9);
    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");
  iw=naxes[1]*9/LINEPIY;
  if(iw==1) wei[0]=1; 

  for(j=0;j<naxes[0];j++){
    buffer[j]=0;
    for(i=0;i<iw;i++) buffer[j]+=z[i+(std-1)*iw][j]*wei[i]/wt[std-1];
    //fprintf(stderr,"%d %g\n",j,buffer[j]);if(j%10==0)getchar();
  }

  naxes[0]=fabs((wlmax-wlmin)/dw)+0.5;
  if(naxes[0]>LINEPIX) naxes[0]=LINEPIX;

  for(i=0;i<naxes[1];i++){
    jw=1;w2=0;
    for(j=0;j<naxes[0];j++){
      zo[i][j]=0;
      no[i][j]=0;
      ws=wlmin+dw*j-dw*0.5;
      we=wlmin+dw*j+dw*0.5;
      while(jw<naxes0-1){
	w1=(w[i/iw][jw-1]+w[i/iw][jw])/2;
	w2=(w[i/iw][jw]+w[i/iw][jw+1])/2;
	if(raw) fac=(w2-w1)/dw;
	else fac=w[i/iw][jw]/expmin*8.36/EFFICIENCY;
	//fprintf(stderr,"%d %d %d %g %g %g %g %g ",i,j,jw,ws,we-ws,w1,w2-w1,zo[i][j]);
	if(w1>we) break;
	else if(w2<ws){
	  //fprintf(stderr,"%g %g\n",z[i][jw],zo[i][j]);
	  jw++;
	  continue;
	}
	else if(w2>we){
	  if(w1>ws){
	    zo[i][j]+=z[i][jw]*(we-w1)/(w2-w1)*fac;
	    no[i][j]+=n[i][jw]*(we-w1)/(w2-w1)*fac*fac;
	  }else{
	  zo[i][j]+=z[i][jw]*(we-ws)/(w2-w1)*fac;
	  no[i][j]+=n[i][jw]*(we-ws)/(w2-w1)*fac*fac;
	  //fprintf(stderr,"%g %g\n",z[i][jw],zo[i][j]);
	  //fprintf(stderr,"%g %g\n",n[i][jw],no[i][j]);
	  }
	  break;
	}
	else if(w1<ws){
	  zo[i][j]+=z[i][jw]*(w2-ws)/(w2-w1)*fac;
	  no[i][j]+=n[i][jw]*(w2-ws)/(w2-w1)*fac*fac;
	}else{ 
	  zo[i][j]+=z[i][jw]*fac;
	  no[i][j]+=n[i][jw]*fac*fac;
	  //fprintf(stderr,"%g %g\n",z[i][jw],zo[i][j]);
	}
	//if(j==831 && i==982)fprintf(stderr,"fac1=%g noise=%g\n",fac,no[i][j]);
	//fprintf(stderr,"%g %g\n",z[i][jw],zo[i][j]);
	//fprintf(stderr,"%g %g\n",n[i][jw],no[i][j]);getchar();
	jw++;
      }
    }
  }

  if(raw==0){
  t1=-99.;t2=4;
  if(argc==8+ot){
    if(strlen(argv[7+ot])==strspn(argv[7+ot],".0123456789")){
      am=atof(argv[7+ot]);
      if(am<1.){
	fprintf(stderr,"Airmass must be > 1. [%lf]\n",am);
	//exit(1);
      }
    }else if(strcmp(argv[7+ot],"NONE")==0 || strcmp(argv[7+ot],"none")==0){
      t1=-1;
      for(j=0;j<LINEPIX;j++) ab[j]=1.;
    }else{
      t1=atof(argv[7+ot]+1);
      if(t1<0 || 10<=t1){
	fprintf(stderr,"Spectral Type must be between 0-9.9 [%f]\n",t1);
	exit(1);
      }
      if(strncmp(argv[7+ot],"F",1)==0) t1+=0;
      else if(strncmp(argv[7+ot],"G",1)==0) t1+=10;
      else if(strncmp(argv[7+ot],"K",1)==0) t1+=20;
      else{
	fprintf(stderr,"Spectral Type must be F, G, or K\n");
	exit(1);
      }
      if(strspn(argv[7+ot]+strlen(argv[7+ot])-1,"0123456789")==1) t2=4;
      else if(strcmp(argv[7+ot]+strlen(argv[7+ot])-2,"IV")!=0 && strcmp(argv[7+ot]+strlen(argv[7+ot])-1,"V")==0) t2=5;
      else if(strcmp(argv[7+ot]+strlen(argv[7+ot])-3,"III")==0) t2=3;
      else{
	fprintf(stderr,"Spectral Type must be III or V or omitted.\n");
	exit(1);
      }
    }
  }
  for(j=0;j<LINEPIX;j++) atm[j]=0.;
  if((fpd=fopen("../atmos.dat","r"))==NULL){
    fprintf(stderr,"Can not open file: ../atmos.dat\n");
    exit(1);
  }
  fprintf(stderr,"Open atmospheric transmission file: ../atmos.dat\n");
  j=0;
  while(!feof(fpd)){
    if(fgets(buf,sizeof(buf),fpd)==NULL) break;
    if(sscanf(buf,"%f %f",&ww,&atmt)!=2){
      fprintf(stderr,"Wrong format\n");
      exit(1);
    }
    if(mode==7){
      for(nj=0;nj<5;nj++){
	jw=(j+2)/4+nj-2;
	if(jw<0) jw=0;
	if(jw>=naxes0) jw=naxes0-1;
	atm[jw]+=atmt*bx[nj]/4;
      }
    }else{
      if(ww<wlmin+dw*j-dw/2) continue;
      if(ww>wlmax-dw/2) break;
      if(ww>wlmin+dw*j+dw/2){
	fprintf(stderr,"No data at %g [um]\n",wlmin+dw*j);
	exit(1);
      }else{
	for(nj=0;nj<5;nj++){
	  jw=j+nj-2;
	  if(jw<0) jw=j;
	  else if(jw>=naxes0) jw=j;
	  atm[jw]+=atmt*bx[nj];
	}
      }
    }
    j++;
  }
  fclose(fpd);
  for(j=0;j<naxes[0];j++){
    ww=wlmin+dw*j;
    //fprintf(stderr,"%g\t%g\t",ww,atm[j]);
    atm[j]=pow(atm[j],am);
    //fprintf(stderr,"%g\n",atm[j]);
  }

  for(i=0;i<naxes1;i++){
    jw=1;w2=0;
    for(j=0;j<naxes[0];j++){
      tr[i][j]=0;
      ws=wlmin+dw*j-dw*0.5;
      we=wlmin+dw*j+dw*0.5;
      while(jw<naxes0-1){
	w1=(w[i][jw-1]+w[i][jw])/2;
	w2=(w[i][jw]+w[i][jw+1])/2;
	if(w1>we) break;
	else if(w2<ws){
	  jw++;
	  continue;
	}
	else if(w2>we){
	  if(w1>ws) tr[i][j]+=buffer[jw]*(we-w1)/(w2-w1)*w[i][jw];
	  else tr[i][j]+=buffer[jw]*(we-ws)/(w2-w1)*w[i][jw];
	  break;
	}
	else if(w1<ws) tr[i][j]+=buffer[jw]*(w2-ws)/(w2-w1)*w[i][jw];
	else tr[i][j]+=buffer[jw]*w[i][jw];
	jw++;
      }
    }
  }

  if(irs==2) {
    ff0[mode-1]*=rf0[mode-1];
    ff1[mode-1]*=rf1[mode-1];
  }
  f0=0;f1=0;f2=0;
  for(j=0;j<naxes[0];j++){
    ww=wlmin+dw*j;
    if(w0min[mode-1]-1e-4<ww && ww<w0max[mode-1]-1e-4) f0+=tr[std-1][j]/atm[j]/((w0max[mode-1]-w0min[mode-1])/dw)*ff0[mode-1];
    if(w1min[mode-1]-1e-4<ww && ww<w1max[mode-1]-1e-4) f1+=tr[std-1][j]/atm[j]/((w1max[mode-1]-w1min[mode-1])/dw)*ff1[mode-1];
    if(w2min[mode-1]-1e-4<ww && ww<w2max[mode-1]-1e-4) f2+=tr[std-1][j]/atm[j]/((w2max[mode-1]-w2min[mode-1])/dw);
  }
  a=(f2-f1)/((w2max[mode-1]+w2min[mode-1])/2-(w1max[mode-1]+w1min[mode-1])/2);
  af0=a/f0;
  fprintf(stderr,"%g %g %g Slope=%g\n",f1,f2,f0,af0);

  if(t2==3) strcpy(fname,"../afIII.dat");
  else if(t2==5) strcpy(fname,"../afV.dat");
  else strcpy(fname,"../afx.dat");
  if((fpd=fopen(fname,"r"))==NULL){
    fprintf(stderr,"Can not open file: %s\n",fname);
      exit(1);
  }
  fprintf(stderr,"Open type_slope file: %s\n",fname);
  for(j=0;j<30;j++){
    for(i=0;i<mode+1;i++) fscanf(fpd,"%s",dmy);
    fscanf(fpd,"%f",&af[j]);
    //fprintf(stderr,"%f\n",af[j]);
    for(i=0;i<7-mode;i++) fscanf(fpd,"%s",dmy);
  }
  fclose(fpd);
    
  if(t1<-90.){
    if(af[29] < af0){
      fprintf(stderr,"##### Selected star is TOO RED ! (Slope > %f)####\n",af[29]);
      daf=af0-af[29];t1=29.;
      fprintf(stderr,"#####   Using K9III + extra Slope %g ... ####\n",daf);
    }
    else if(af0 < af[0]){
      fprintf(stderr,"##### Selected star is TOO BLUE ! (Slope < %f)####\n",af[0]);
      daf=af0-af[0];t1=0.;
      fprintf(stderr,"#####   Using F0V + extra Slope %g ... ####\n",daf);
    }
    else{
      for(j=0;j<29;j++){
	if(af[j]<=af0 && af0<af[j+1]) t1=(j*(af[j+1]-af0)+(j+1)*(af0-af[j]))/(af[j+1]-af[j]);
      }
      fprintf(stderr,"Spectral Type: ");
      if(0<=t1 && t1<10) fprintf(stderr,"F%3.1f",t1);
      else if(10<=t1 && t1<20) fprintf(stderr,"G%3.1f",t1-10);
      else if(20<=t1 && t1<30) fprintf(stderr,"K%3.1f",t1-20);
      else fprintf(stderr,"---\n");
      if(0<=t1 && t1<=15) fprintf(stderr,"V\n");
      else if(15<t1 && t1<21) fprintf(stderr,"V-III\n");
      else if(21<=t1 && t1<30) fprintf(stderr,"III\n");
    }
  }else{
    afm=af[(int)t1]+(t1-(int)t1)*(af[(int)t1+1]-af[(int)t1]);
    ff1t=(f2-afm*f0*((w2max[mode-1]+w2min[mode-1])/2-(w1max[mode-1]+w1min[mode-1])/2))/f1;
    fprintf(stderr,"Model f1, Slope=%f, %f\n",ff1t*f1,afm);
    if(irs==2) {ff1[mode-1]/=rf1[mode-1];rf1[mode-1]*=ff1t;}
    else ff1[mode-1]*=ff1t;
    if((fpd=fopen("conv_w_fnu.dat","w"))!=NULL){
      fprintf(stderr,"Save color correction file: conv_w_fnu.dat\n");
      for(i=0;i<7;i++) fprintf(fpd," %f",ff1[i]);
      fprintf(fpd,"\n");
      for(i=0;i<7;i++) fprintf(fpd," %f",rf1[i]);
      fprintf(fpd,"\n");
      fclose(fpd);
    }
  }

  if(-0.1<t1 && t1<29.1){
    if(t2==3){
      if(t1<15) strcpy(fname,"../irtflib_FIII.dat");
      else if(t1<25) strcpy(fname,"../irtflib_K1III.dat");
      else strcpy(fname,"../irtflib_K4III.dat");
      if(t1<5) strcpy(fname2,"../irtflib_FIII.dat");
      else if(t1<21) strcpy(fname2,"../irtflib_GIII.dat");
      else strcpy(fname2,"../irtflib_K4III.dat");
    }else if(t2==5){
      if(t1<15) strcpy(fname,"../irtflib_FV.dat");
      else strcpy(fname,"../irtflib_KV.dat");
      if(t1<5) strcpy(fname2,"../irtflib_FV.dat");
      else if(t1<23.5) strcpy(fname2,"../irtflib_GV.dat");
      else strcpy(fname2,"../irtflib_KV.dat");
    }else{
      if(t1<15) strcpy(fname,"../irtflib_FV.dat");
      else if(t1<25) strcpy(fname,"../irtflib_K1III.dat");
      else strcpy(fname,"../irtflib_K4III.dat");
      if(t1<5) strcpy(fname2,"../irtflib_FV.dat");
      else if(t1<21) strcpy(fname2,"../irtflib_GV.dat");
      else strcpy(fname2,"../irtflib_K4III.dat");
    }
    for(j=0;j<LINEPIX;j++) ab[j]=0.;
    if((fpd=fopen(fname,"r"))==NULL){
      fprintf(stderr,"Can not open file: %s\n",fname);
      exit(1);
    }
    fprintf(stderr,"Open stellar absorption file: %s\n",fname);
    if((fps=fopen(fname2,"r"))==NULL){
      fprintf(stderr,"Can not open file: %s\n",fname2);
      exit(1);
    }
    fprintf(stderr,"Open stellar absorption file: %s\n",fname2);
    j=0;
    while(!feof(fpd)&&!feof(fps)){
      if(fgets(buf,sizeof(buf),fpd)==NULL) break;
      if(sscanf(buf,"%f %f",&ww,&abt)!=2){
	fprintf(stderr,"Wrong format\n");
	exit(1);
      }
      if(fgets(buf,sizeof(buf),fps)==NULL) break;
      if(sscanf(buf,"%f %f",&ww,&abt2)!=2){
	fprintf(stderr,"Wrong format\n");
	exit(1);
      }
      if(t2==5){
	if(5<=t1 && t1<15) abt=(abt*(15-t1)+abt2*(t1-5))/10;
	else if(15<=t1 && t1<23.5) abt=(abt2*(23.5-t1)+abt*(t1-15))/8.5;
      }else{
	if(5<=t1 && t1<15) abt=(abt*(15-t1)+abt2*(t1-5))/10;
	else if(15<=t1 && t1<21) abt=(abt2*(21-t1)+abt*(t1-15))/6;
	else if(21<=t1 && t1<25) abt=(abt*(25-t1)+abt2*(t1-21))/4;
      }

      if(mode==7){
	for(nj=0;nj<5;nj++){
	  jw=(j+2)/4+nj-2;
	  if(jw<0) jw=0;
	  if(jw>=naxes[0]) jw=naxes[0]-1;
	  ab[jw]+=abt*bx[nj]/4;
	}
      }else{
	if(ww<wlmin+dw*j-dw/2) continue;
	if(ww>wlmax-dw/2) break;
	if(ww>wlmin+dw*j+dw/2){
	  fprintf(stderr,"No data at %g [um]\n",wlmin+dw*j);
	  exit(1);
	}else{
	  for(nj=0;nj<5;nj++){
	    jw=j+nj-2;
	    if(jw<0) jw=j;
	    else if(jw>=naxes[0]) jw=j;
	    ab[jw]+=abt*bx[nj];
	  }
	}
      }
      j++;
    }
    fclose(fpd);
    fclose(fps);
    af3=1.228519925e-03*t1*t1+2.651552416e-02*t1-6.757720080e-01;
    af5=3.598987915e-04*t1*t1+2.887846009e-02*t1-6.714439294e-01;
    if(t2==3) a=af3*f0;
    else if(t2==5) a=af5*f0;
    else{
      if(t1<=15) a=(af5+daf)*f0;
      else if(21<=t1) a=(af3+daf)*f0;
      else a=(af5*(21-t1)+af3*(t1-15))/6*f0;
    }
  }
  fprintf(stderr,"Intrinsic Slope=%g\n",a/f0);
  
  for(j=0;j<naxes[0];j++){
    ww=wlmin+dw*j;
    for(i=0;i<naxes1;i++) tr[i][j]*=(wt[i]/(a*(ww-(w0max[mode-1]+w0min[mode-1])/2)+f0)/ab[j]);
    //fprintf(stderr,"%d %g\n",j,tr[100][j]);if(j%10==0)getchar();
  }

  if(mode==7){
    trs1=0.;trs2=0.;
    for(j=100;j<300;j++) trs1+=tr[100][j];
    for(j=700;j<900;j++) trs2+=tr[100][j];

    fprintf(stderr,"Efficiency ratio (1.0um/1.3um) = %g (Typical value: 0.094)\n",trs1/trs2);
    if(trs1/trs2*10.> 1.3){
      fprintf(stderr,"##### Reddening correction is required. ####\n");
    }
  }

  for(i=0;i<naxes[1];i++){
    for(j=0;j<naxes[0];j++){
      fac1=0.;fac2=0.;fac3=0.;
      for(jw=j-2;jw<=j+2;jw++){
	if(jw<0) fac=tr[i/iw][0];
	if(jw>=naxes[0]) fac=tr[i/iw][naxes[0]-1];
	else fac=tr[i/iw][jw];
	if(fac>fac1){fac3=fac2;fac2=fac1;fac1=fac;}
	else if(fac>fac2){fac3=fac2;fac2=fac;}
	else if(fac>fac3){fac3=fac;}
	//fprintf(stderr,"%d %d %g %g %g %g\n",j,jw,fac,fac1,fac2,fac3);getchar();
      }
      fac=fac3;
      //fac=tr[i/iw][j];
      if(fac<1e-5) fac=1e-5;
      zo[i][j]*=irs/fac;
      no[i][j]*=irs*irs/fac/fac;
      //if(j==831 && i==982)fprintf(stderr,"fac2=%g irs=%d noise=%g\n",fac,irs,no[i][j]);
    }
  }
  }

  fprintf(stderr,"Write %s  ... ",argv[4]);
  status=0;
  if (fits_create_file(&fp, argv[4], &status))
    printerror(status);
  if (fits_create_img(fp, bitpix, naxis, naxes, &status))
    printerror(status);
  if (fits_write_key(fp, TSTRING, "CTYPE1", "WAVELENGTH", "", &status))
    printerror(status);
  if (fits_write_key(fp, TSTRING, "CUNIT1", "Angstrom", "", &status))
    printerror(status);
  if (fits_write_key(fp, TSTRING, "CTYPE2", "SPECTRUM", "", &status))
    printerror(status);
  if (fits_write_key(fp, TSTRING, "CUNIT2", "Number", "", &status))
    printerror(status);
  tmpval=1.;
  if (fits_write_key(fp, TFLOAT, "CRPIX1", &tmpval, "", &status))
    printerror(status);
  tmpval=wlmin*10000;
  if (fits_write_key(fp, TFLOAT, "CRVAL1", &tmpval, "", &status))
    printerror(status);
  tmpval=dw*10000;
  if (fits_write_key(fp, TFLOAT, "CD1_1", &tmpval, "", &status))
    printerror(status);
  tmpval=5.;
  if (fits_write_key(fp, TFLOAT, "CRPIX2", &tmpval, "", &status))
    printerror(status);
  tmpval=1.;
  if (fits_write_key(fp, TFLOAT, "CRVAL2", &tmpval, "", &status))
    printerror(status);
  tmpval=1./iw;
  if (fits_write_key(fp, TFLOAT, "CD2_2", &tmpval, "", &status))
    printerror(status);
  fpixel=1;
  nelements=naxes[0];
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  for(i=0;i<naxes[1];i++){
    if (fits_write_img(fp, TFLOAT, fpixel, nelements, zo[i], &status))
      printerror(status);
    fpixel+=naxes[0];
  }
  if (fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  if(raw==0){
  sscanf(argv[4],"%s",fname);
  fname[strlen(argv[4])-5]='\0';
  sprintf(fname,"%s_nosq.fits",fname);
  remove(fname);
  fprintf(stderr,"Write %s  ... ",fname);
  status=0;
  if (fits_create_file(&fp, fname, &status))
    printerror(status);
  if (fits_create_img(fp, bitpix, naxis, naxes, &status))
    printerror(status);
  fpixel=1;
  nelements=naxes[0];
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  for(i=0;i<naxes[1];i++){
    if (fits_write_img(fp, TFLOAT, fpixel, nelements, no[i], &status))
      printerror(status);
    fpixel+=naxes[0];
  }
  if (fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  naxes[1]=naxes1;
  fprintf(stderr,"Write check.fits  ... ");
  status=0;
  if (fits_create_file(&fp, "check.fits", &status))
    printerror(status);
  if (fits_create_img(fp, bitpix, naxis, naxes, &status))
    printerror(status);
  fpixel=1;
  nelements=naxes[0];
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  for(i=0;i<naxes[1];i++){
    if (fits_write_img(fp, TFLOAT, fpixel, nelements, tr[i], &status))
      printerror(status);
    fpixel+=naxes[0];
  }
  if (fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\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 */
}