#include <stdio.h>
#define nbrowte 4160*4096*4
char head[72][80];
char a[nbrowte];
char f1[30],f2[30];
FILE *fp;	
int ibit,n1,n2,ipos;
int bcol,lcol,brow,lrow;
double x;
main(ac,av)
int ac; char *av[];                        
{
  int i;			
  if(ac<7){
    printf("\tUsage:  substract infile outputfile bcol lcol brow lrow\n\n");
    printf("\t        begin_col, length_col,  begin_row, length_row\n");
    printf("\tEx:     substract p*.fit p2048.fit 1024 2048 1024 2048\n");
    printf("\t        p2048.fit is a 2k*2k file of origin center\n");
    printf("\t\t\t\t\t\t\t  60321\n");
    exit(0);
  }  
  sscanf(&av[3][0],"%d",&bcol);  sscanf(&av[4][0],"%d",&lcol);
  sscanf(&av[5][0],"%d",&brow);  sscanf(&av[6][0],"%d",&lrow);
  fp=fopen(av[1],"rb"); if(fp==0){ printf("\n\tinfile not found!\n"); exit(0); }
  fread(head,80,72,fp);
  ipos=indexpos(head,"BITPIX  ",72);
  if(ipos==72){ printf("\n\tinfile not a fits file!\n"); exit(0); }
  sscanf(&head[ipos][23],"%d",&ibit);	  
  ipos=indexpos(head,"NAXIS1  ",72);
  sscanf(&head[ipos][23],"%d",&n1);	  
  sscanf(&head[ipos+1][23],"%d",&n2);	  
  if(ibit<0)ibit=-ibit; ibit/=8;
  printf("Size: %d %d %d\n",ibit,n1,n2);
  if(bcol+lcol>n1){ printf("bcol+lcol>n1!\n"); exit(0); }
  if(brow+lrow>n2){ printf("brow+lrow>n2!\n"); exit(0); }
//  printf("%s %s %d %d %d %d %d %d\n",f1,f2,n1,n2,bcol,lcol,brow,lrow);
  fread(a,n1*n2,ibit,fp); fclose(fp);
  if(ibit==2)ch2(a,a,n1,n2,bcol,lcol,brow,lrow);
        else ch4(a,a,n1,n2,bcol,lcol,brow,lrow);
  sprintf(&head[ipos][22],"%8d",lcol); head[ipos][30]=32;
  sprintf(&head[ipos+1][22],"%8d",lrow); head[ipos+1][30]=32;
  ipos=indexpos(head,"CRPIX1  ",72);
  sscanf(&head[ipos][10],"%le",&x);
  x-=bcol;
  putx();
  ipos=indexpos(head,"CRPIX2  ",72);
  sscanf(&head[ipos][10],"%le",&x);
  x-=brow;
  putx();
	
  fp=fopen(av[2],"wb");
  fwrite(head,80,72,fp);
  fwrite(a,lcol*lrow,ibit,fp);
  fclose(fp);
  printf("\nok!\n");
}

ch4(a,b,n1,n2,bcol,lcol,brow,lrow)
float a[n2][n1],*b;
int n1,n2,bcol,lcol,brow,lrow;
{
  int j,i;
  for(j=0;j<lrow;j++)for(i=0;i<lcol;i++)*b++=a[j+brow][i+bcol];
}  

ch2(a,b,n1,n2,bcol,lcol,brow,lrow)
short a[n2][n1],*b;
int n1,n2,bcol,lcol,brow,lrow;
{
  int j,i;
  for(j=0;j<lrow;j++)for(i=0;i<lcol;i++)*b++=a[j+brow][i+bcol];
}  

putx()
{
  char c[20];
  int i;
  for(i=0;i<20;i++)c[i]=32;
  sprintf(c,"%le",x);
  for(i=0;i<20;i++)if(c[i]==0)c[i]=32;
  for(i=0;i<20;i++)head[ipos][10+i]=c[i];
}