/* gcc -o readfits readfits.c  */
#include <stdio.h>
#define maxbyte 2048*2048*4    /* float arrays 2048,2048 */

float bscale,bzero0;
char  head[36][80];
char  image[maxbyte];

main(ac,av)
int ac; char *av[];
{
  FILE  *fp;
  int   i,n,ipos,n1,n2,bitpix;
 
  if(ac<2){ printf("\n\t Usage: readfits fitsfile \n\n"); exit(0); }
  fp=fopen(av[1],"rb");
  if(fp==0){ printf("\n\t fitsfile not found!\n"); exit(0); }
  
/* read fits head */
  fread(head,1,2880,fp);
  if(indexpos(head,"SIMPLE  ",36)!=0)goto error;
  if(indexpos(head,"BITPIX  ",36)!=1)goto error;
  sscanf(&head[1][26],"%d",&bitpix);
  if((ipos=indexpos(head,"NAXIS1  ",36))==36)goto error;
  sscanf(&head[ipos][26],"%d",&n1);
  if((ipos=indexpos(head,"NAXIS2  ",36))==36)goto error;
  sscanf(&head[ipos][26],"%d",&n2);
  bscale=1.; bzero0=0.; i=0;
readhead:
  i++;  if(i>5)printf("\n I cannot finf keywords END (with 5 space)\n\n");
        if(i>5)exit(0);
  if((ipos=indexpos(head,"BSCALE  ",36))!=36)
    sscanf(&head[ipos][20],"%f",&bscale);
  if((ipos=indexpos(head,"BZERO   ",36))!=36)
    sscanf(&head[ipos][20],"%f",&bzero0);
  if(indexpos(head,"END     ",36)==36){
    fread(head,1,2880,fp);
    goto readhead;
  }

  printf("datatype: ");
  if(bitpix==16)printf(" short image[%d][%d]\n",n2,n1);
  if(bitpix==32)printf("  long image[%d][%d]\n",n2,n1);
  if(bitpix==-32)printf("float image[%d][%d]\n",n2,n1);

/* read fits body */
  n=n1*n2*4; if(bitpix==16)n/=2;
  fread(image,1,n,fp);
  fclose(fp);
  i=0; if(bscale==1. && bzero0==0. )i=1;
  if(bitpix== 16){ swap2(image,n); if(i==0)chshort(image,n1,n2); } 
  if(bitpix== 32){ swap4(image,n); if(i==0) chlong(image,n1,n2); }
  if(bitpix==-32){ swap4(image,n); if(i==0)chfloat(image,n1,n2); }
    
/* you can use data as above subroutines according data_type */

  printf("ok !\n"); 
  exit(0);

error:
  printf("\n\t Not a 2D fits file!\n"); exit(0);
}

chshort(a,n1,n2)
int n1,n2;
short a[][n1];
{
  int i,j;
  for(j=0;j<n2;j++)for(i=0;i<n1;i++)a[j][i]=bscale*a[j][i]+bzero0;
}
  
chlong(a,n1,n2)
int n1,n2;
long a[][n1];
{
  int i,j;
  for(j=0;j<n2;j++)for(i=0;i<n1;i++)a[j][i]=bscale*a[j][i]+bzero0;
}
  
chfloat(a,n1,n2)
int n1,n2;
float a[][n1];
{
  int i,j;
  for(j=0;j<n2;j++)for(i=0;i<n1;i++)a[j][i]=bscale*a[j][i]+bzero0;
}
  
indexpos(head,f1,n)
char head[][80],f1[8];
int n;
{
  int i,j;
  for(i=0;i<n;i++){
    for(j=0;j<8;j++) if(head[i][j]!=f1[j])goto l10;
    return i;
l10:
    continue;
  }
  return i;
}

swap2(char *a, int n)
{
  register char *p;
  register int i;
  register char c;
  p = a;
  for(i=0;i<n;i+=2){
    c = *p;
    *p++ = *(p+1);
    *p++ = c;
  }
}

swap4(char *a, int n)       /*     long a[n] or float a[n]    */
{
  register char *p;
  register int i;
  register char c,d;
  p = a;
  for(i=0;i<n;i+=4){
    c = *p;
    *p++ = *(p+3);
    d = *p;
    *p++ = *(p+1);
    *p++ = d;
    *p++ = c;
  }
}