#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

void jiang_()
{
 printf("\n********** C subroutine ***************************\n");
 printf("long file_size( name )             char *name\n");
 printf("int  index( a,c )                  char a[],c\n");
 printf("int  nindex( a,b )                 char a[],b[]\n");
 printf("void stop( s )                     char s[]\n");
 printf("int  indexpos( head,f1,n )         char head[][80],f1[8]\n");
 printf("                                   int n; (input, 36 or 72..)\n");
 printf("int  c_indexpos(head,f1);          return 0--72\n");
 printf("void swap4( a,nbyte )              float a[] or long a[]; int n\n");
 printf("void swap2( a,nbyte )              short a[]; int n\n");
 printf("void white_black(map,n1,n2,peak,sigm)\n");
 printf("void histat(ihist,xpeak,sigma,n)\n");
 printf("void ad_lb(al4,de4,pl,pb)\n");
 printf("void toms(aa,cc,k)                 float aa; char cc[12]; k=1,0,2\n");
 printf("void usrid_(a)                     char a[]\n");
 printf("void fdate_(a)                     char a[]\n");
 printf("void filesize_(name,flen)\n");
 printf("void tvcoordc_(f1,a)               write a[72][80] to fits_head\n");
 printf("void tvread_(f1,a,len,pos)\n");
 printf("void tvwrite_(f1,a,len,pos)\n");
 printf("void dispdot_()\n");
 printf("void readfits_(f1,a,n1,n2)         char f1[],a[], *n1,*n2\n");
 printf("void readc_(f1,map,ix,iy)\n");
 printf("\n************* Fortran routine ***********\n");
 printf("pgxgray(map,n1,n2,nx1,nx2,ny1,ny2,white,black,irat)   irat=1,2\n");
 printf("pgximage(map,n1,n2,white,black)\n");
 printf("pgximage0(map,n1,n2,white,black)\n");
 printf("pgxchinese(x,y,ch,n,j1,j2)             i4_ch(n), j2 for direct\n");
 printf("turnx2(a,n1,n2,it)              i*2_a(n1,n2);it=90,...,-360,..\n");
 printf("readxfile(file,ihead,map,nbyte)                   b_map(nbyte)\n");
 printf("whitexblack(map,n1,n2,peak,sigma)                i2_map(n1,n2)\n");
 printf("histat(i4_ihist,mode,maxh,xmean,xpeak,fsigma,gsigma,n)\n");
 printf("parbol(a,n,xpeak,sigma)                    least 2, up to x**4\n");
 printf("toms(a,c,k)              r_a, c*12_a, k=1 for hour; =0 for deg\n");
 printf("itohd(c,alpha,delta,epoch)           c*(*)_c, is a readin line\n");
 printf("astprs(r1,d1,e1,r2,d2,e2)                 1 epoch to 2, in h,d\n");
 printf("standc(rc,dc,ra,de,xi,xn)                         in arc, r->x\n");
 printf("astand(rc,dc,xi,xn,ra,de)                         in arc, x->r\n");
 printf("xytoad(x,a)                      2_6 convert, r8_x(2,6),a(2,6)\n");
 printf("plate(ad,xy,w,n,c,nc)             r8_ad,xy(2,n)),c(2,6) b_w(n)\n");
 printf("slove or solve8(a,b,m)                    a(6,6),b(6) b-in-out\n");
 printf("correct216(xi,xn)                                       in arc\n");
 printf("swap2(a,nbyte)   swap4(a,nbyte)\n");
 printf("lnblnk(string)\n");
 printf("unlink(file)\n");
 printf("iand(i2,j2), ior(i2,j2), inot(i2)\n");
}

long file_size( name )
char *name;
{
     long eof_ftell;
     FILE *file;

     file = fopen( name, "r");
     if ( file == NULL ) return( 0L );
     fseek( file, 0L, 2);
     eof_ftell = ftell( file );
     fclose( file );
     return( eof_ftell );
}

index(a,c)
char a[],c;
{
  int i,j;
  j=strlen(a);
  for(i=0;i<j;i++)if(a[i]==c)break;
  return (i<j)?i:-1;
}
 
nindex(a,b)
char a[],b[];
{
  int i,j,k;
  char c;
  j=strlen(a);
  k=strlen(b);
  c=b[0];
  for(i=0;i<j;i++)if(a[i]==c && strncmp(&a[i],b,k)==0)break;
  return (i<j)?i:-1;
}
 
stop(s)
char s[];
{
  printf("%s\n",s); exit(0);
}
 
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(a,n)
char a[]; int n;
{
  int i;
  char c;
  for(i=0;i<n;i+=2){
    c=a[i];
    a[i]=a[i+1];
    a[i+1]=c;
  }
}

swap4(a,n)
char a[]; int n;
{
  int i;
  char c;
  for(i=0;i<n;i+=4){
    c=a[i];
    a[i]=a[i+3];
    a[i+3]=c;
    c=a[i+1];
    a[i+1]=a[i+2];
    a[i+2]=c;
  }
}
*/

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;
  }
}

white_black(map,n1,n2,peak,sigma)
short map[];
int n1,n2;
float *peak,*sigma;
{
  int *ihist;
  ihist=(int*)malloc(32767*4);
  white_black_1(map,n1,n2,peak,sigma,ihist);
  free(ihist);
}

white_black_1(map,n1,n2,xpeak,sigma,ihist)
short map[]; int n1,n2;
float *xpeak, *sigma;
long  ihist[];
{
  int i,ii;
  int imax=0;
  for(i=0;i<32767;i++)ihist[i]=0;
  for(i=0;i<n1*n2;i++){
    ii=map[i];
    ii+=101;
    if(ii<0 || ii >= 32767)continue;
    ihist[ii]++;
    if(ii>imax)imax=ii;
  }
  histat(ihist,xpeak,sigma,imax);
  *xpeak-=100.0;
}

histat(ihist,xpeak,fsigma,n)
long ihist[];
float *xpeak,*fsigma;
int   n;
{
  int ihtot=0, ihsum=0, iwid=0, maxh=0, mode=0;
  int i,ii,j,jj,k,m,nn, isigm, nfilt, noff, ixpeak, ilow, ihih,ilim,icount;
  float xpar[101],buf[201];
  float xmean=0., xmed=0, sigmsq=0., xcount=0;
  float rttpi, conv, xx, temp, xp,fs,gs,cogd,cogn,cog,xnumb,sd;

  rttpi=2.506628275;
  for(i=1;i<=n;i++){
    nn=n+1-i;
    if(ihist[nn]!=0)break;
  }
  for(i=1;i<=nn;i++)ihtot+=ihist[i];
  for(i=1;i<=nn;i++){
    ihsum+=ihist[i];
    if(ihsum < 0.1*ihtot || ihsum > 0.9*ihtot)continue;
    if(ihist[i] <= maxh)continue;
    maxh=ihist[i];
    mode=i;
  }
  for(i=1;i<=nn;i++){
    temp=ihist[i];
    xcount+=temp;
    xmean+=temp*i;
    sigmsq+=temp*i*i;
  }
  for(i=1;i<=nn;i++){
    xmed+=ihist[i];
   if(xmed > 0.5*xcount){ xmed=i; break; }
  }
  xmean/=xcount;
  gs=xcount/(maxh*rttpi);
  if(gs > 0.1*nn)gs=0.1*nn;
  isigm=gs+0.5;
  if(isigm<3)isigm=3;
  if(isigm>50)isigm=50;
  nfilt=isigm*4+1;
  noff=(nfilt+1)/2;
  conv=0.5/(float)(isigm*isigm);
  for(i=1;i<=noff;i++){
    xx=i-noff;
    buf[i]=exp(-conv*xx*xx);
    buf[nfilt+1-i]=buf[i];
  }
  xp=mode;
  for(k=1;k<=2;k++){
    ixpeak=xp+0.5;
    if(ixpeak > nn)ixpeak=nn;
    if(ixpeak <1 )ixpeak=1;
    ilow=ixpeak-isigm;
    if(ilow <1 )ilow=1;
    ihih=ixpeak+isigm;
    if(ihih >nn)ihih=nn;
    m=ihih-ilow+1;
    ii=1;  cogd=cogn=0.;
    for(i=ilow;i<=ihih;i++){
      cogd+=ihist[i];
      cogn+=ihist[i]*i;
      temp=xnumb=0.;
      for(j=1;j<=nfilt;j++){
	jj=j-noff;
	if(i+jj < 1 || i+jj > nn)continue;
	temp+=ihist[i+jj]*buf[j];
	xnumb+=buf[j];
      }
      temp = (xnumb > 0.) ? temp/=xnumb : 0. ;
      if(temp < 1.)temp=1.;
      xpar[ii++]=log((double)temp);
    }
    parbol(xpar,m,&xp,&sd);
    xx=sd*sd-isigm*isigm;
    if(xx < 0.)xx=0.;
    fs=sqrt((double)xx);
    xp+=ilow-1;
    cog=(cogd > .5) ? cogn/cogd : xp;
    xx=xp-mode; if(xx <0.)xx=-xx;
    if(xx > gs)xp=cog;
  }
  mode--;
  xp=xp-1.0;
  xmean=xmean-1.0;
  xmed=xmed-1.0;
  if(xp > nn-gs)xp=xmed;
  if(xp < gs)xp=xmed;
  ixpeak=xp+1.5;
  icount=ihist[ixpeak];
  ilim=0.68*xcount+0.5;
  for(i=1;i<=nn;i++){
    ii=ixpeak+i;
    if(ii <= nn){
      iwid++;
      icount+=ihist[ii];
      if(icount > ilim)break;
    }
    ii=ixpeak-i;
    if(ii >= 1){
      iwid++;
      icount+=ihist[ii];
      if(icount > ilim)break;
    }
  }
  gs=0.5*iwid;
  if(gs<fs)fs=gs;
  if(fs<1.)fs=1.;
  *xpeak=xp;
  *fsigma=fs;
}

parbol(a,n,xpeak,sd)
float a[],*xpeak,*sd; int n;
{
  float aa,fi,fi2;
  float cmatx[25][25],bvect[25];   /* same as fortran , 1--24 */
  float	suma=0.,sumay=0.,sumayy=0.,sumy1=0.,sumy2=0.,sumy3=0.,sumy4=0.;
  int i,j;

  for(i=1;i<25;i++){
    bvect[i]=0.;
    for(j=1;j<25;j++)cmatx[i][j]=0.;
  }
  for(i=1;i<=n;i++){
    aa=a[i]; fi=i; fi2=i*i;
    suma+=aa;  sumay+=aa*fi; sumayy+=aa*fi2;
    sumy1+=fi; sumy2+=fi2;   sumy3+=fi*fi2;   sumy4+=fi2*fi2;
  }
  cmatx[1][1]=n;     cmatx[1][2]=sumy1; cmatx[1][3]=sumy2;
  cmatx[2][1]=sumy1; cmatx[2][2]=sumy2; cmatx[2][3]=sumy3;
  cmatx[3][1]=sumy2; cmatx[3][2]=sumy3; cmatx[3][3]=sumy4;
  bvect[1]=suma;     bvect[2]=sumay;    bvect[3]=sumayy;
  solve(cmatx,bvect,3);
  if(bvect[3] == 0) *xpeak= *sd=0.;
  else {
    aa=bvect[3]*2.;
    *xpeak=-bvect[2]/aa;
    if(aa < 0.)aa=-aa;
    *sd=sqrt((double)1.0/aa);
  }
}

solve(a,b,m)
float a[25][25],b[25]; int m;
{
  int i,j,k,l;
  float big,temp;

  for(i=1;i<m;i++){
    big=0.;
    for(k=i;k<=m;k++){
      temp=a[k][i];
      if(temp < 0.)temp=-temp;
      if(temp > big){ big=temp; l=k; }
    }  
    if(big == 0.){
      for(k=1;i<=m;k++)b[k]=0;
      printf(" Zero determinant from SOLVE\n");
      return;
    }
    if(i != l){
      for(j=1;j<=m;j++){ temp=a[i][j]; a[i][j]=a[l][j]; a[l][j]=temp;}
      temp=b[i]; b[i]=b[l]; b[l]=temp;
    }
    big=a[i][i];
    for(j=i+1;j<=m;j++){
      temp=a[j][i]/big;
      b[j]-=temp*b[i];
      for(k=i;k<=m;k++)a[j][k]-=temp*a[i][k];
    }
  }
  for(i=1;i<=m;i++){ 
    l=m+1-i;
    if(a[l][l]==0.){ b[l]=0.; continue; }
    temp=b[l];
    if(l != m)for(j=2;j<=i;j++){ k=m+2-j; temp-=a[l][k]*b[k]; }
    b[l]=temp/a[l][l];
  }
}

ad_lb(al4,de4,pl,pb)
float al4,de4,*pl,*pb;
{
  double al,de,cy,xl,xm,xn;
  cy=atan(1.)/45.;
  al=al4; de=de4;
  xl=cos(de)*cos(al);
  xm=0.9175*cos(de)*sin(al)+0.3978*sin(de);
  xn=-.3978*cos(de)*sin(al)+0.9175*sin(de);
  *pl=atan2(xm,xl)/cy;
  if(*pl < 0.)*pl+=360.;
  *pb=asin(xn)/cy;
}

toms(aa,cc,k)
/*
c hour or degree to char_line
c if k=1 (hour case) in **:**:**.*
c    k=0 (degree       -**:**:**
c    k=2 (degree       ***:**:**
*/
char cc[12];
float aa;
int k;
{
  int i,j,m,n;
  float a,x;
  a=aa; if(a < 0.)a=-a;
  i=a; x=(a-i)*60.;
  j=x; x=(x-j)*60.;
  m=x;
  n=(x-m)*10.+0.5;
  if(k == 0)   n=n+5;
  if(n >= 10) { m++;  n=0; }
  if(m == 60) { j++;  m=0; }
  if(j == 60) { i++;  j=0; }
  if(k !=2)sprintf(cc," %02d:%02d:%02d.%d\0",i,j,m,n);
  if(k ==2)sprintf(cc,"%3d:%02d:%02d.%d\0",i,j,m,n);
  if(k !=2)cc[0]=(aa<0.)?'-':' ';
  if(k!=1)cc[9]=0;
}

usrid_(a)
char a[];
{
  char b[26];
  struct tm *t;
  time_t it;
  it=time(NULL);
  t=localtime(&it);
  strcpy(a,getenv("USER"));
  strcpy(b,asctime(t));
  b[0]=b[1]=b[2]='-';
  strcat(a,b);
}

fdate_(a)
char a[];
{
  struct tm *t;
  time_t it;
  it=time(NULL);
  t=localtime(&it);
  strcpy(a,asctime(t));
}

void filesize_( name, flen )
char *name;
long *flen;
{
     FILE *fp;
     char a[80];
     int i;
     strcpy(a,name);
     for(i=0;i<79;i++)if(a[i]==32)break;
     a[i]=0;
         
     fp = fopen( a, "r");
     *flen=0L;
     if ( fp == NULL ) return;
     fseek( fp, 0L, 2);
     *flen = ftell( fp );
     fclose( fp );
     return;
}

void tvcoordc_(f1,a)
char f1[],a[];
{
  FILE *fp;
  if((fp=fopen(f1,"r+b"))==NULL)printf("error\n");
  fseek(fp,0l,0);
  fwrite(a,80,72,fp);
  fclose(fp);
}

void tvread_(f1,a,len,pos)
char f1[],a[];
int *len,*pos;
{
  FILE *fp;
  int i,j;
  i=*len; j=*pos;
  if((fp=fopen(f1,"r+b"))==NULL)printf("error\n");
  fseek(fp,j,0);
  fread(a,1,i,fp);
  fclose(fp);
}

void tvwrite_(f1,a,len,pos)
char f1[],a[];
int *len,*pos;
{
  FILE *fp;
  int i,j;
  i=*len; j=*pos;
  if((fp=fopen(f1,"r+b"))==NULL)printf("error\n");
  fseek(fp,j,0);
  fwrite(a,1,i,fp);
  fclose(fp);
}

void writefits_(f1,a,len,ip)
char f1[],a[];
int *len,*ip;
{
  FILE *fp;
  int i,j;
  if(*ip==1){
    fp=fopen(f1,"wb");
    fwrite(a,1,*len,fp);
    fclose(fp);
  } else {
    fp=fopen(f1,"r+b");
    fseek(fp,0,2);
    fwrite(a,1,*len,fp);
    fclose(fp);
  }
}

void dispdot_()
{
  printf("."); fflush(stdout);
}

c_indexpos(head,f1)
char head[72][80],f1[8];
{
  int i,j;
  for(i=0;i<72;i++){
    for(j=0;j<8;j++) if(head[i][j]!=f1[j])goto l10;
    return i;
l10:
    continue;
  } 	
  return i;
}

  char small[24]={27,'(','s','1','7','h','9','v','0','t',
       27,'&','l','8','D',27,'&','a',57,'l','6','4','R',13};
  char tail[5]={27,'&','l','0','H'};
void writetail_(head)
char head[72][80];
{
  char c; int i,j;
  FILE *fp;
  fp=fopen("pgplot.hp","r+b");
  fseek(fp,-5,2);
  fwrite(small,1,24,fp);
  i=c_indexpos(head,"DATE-OBS");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"TIME    ");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"EXPOSURE");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"RA      ");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"DEC     ");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"HA      ");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"EPOCH   ");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"OBJECT  ");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"OBSERVER");
  if(i<72){ fwrite(&head[i][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  i=c_indexpos(head,"GALLONG ");
  if(i<72)for(j=0;j<11;j++){
    fwrite(&head[i++][0],1,79,fp); fputc(13,fp); fputc(10,fp);}
  fwrite(tail,1,5,fp);
  fclose(fp);
}

void readfits_(f1,a,n1,n2)
char f1[],a[];
int *n1,*n2;
{
  FILE *fp;
  int i,j;
  char b[36][80],f2[60];
  for(i=0;i<60;i++){
    if(f1[i]<=' '){ f2[i]=0; break; }
    f2[i]=f1[i];
  }
  *n1=*n2=0;
  if((fp=fopen(f2,"r+b"))==NULL){ printf("open error\n"); return; }
  fread(b,80,36,fp);
  i=indexpos(b,"BITPIX  ",36); if(i==36)return;
  sscanf(&b[i][23],"%d",&j);   if(j<0)j=-j; j/=8; 
  i=indexpos(b,"NAXIS1  ",36); if(i==36)return;
  sscanf(&b[i+0][23],"%d",n1);
  sscanf(&b[i+1][23],"%d",n2);
  if(indexpos(b,"END     ",36)!=36)goto l11;
  i=0;
l10:
  i++; if(i>9){ *n1=*n2=0; return; }
  fread(b,36,80,fp);
  if(indexpos(b,"END     ",36)==36)goto l10;
l11:
  i=j*(*n1)*(*n2);
  fread(a,1,i,fp);
  fclose(fp);
  if(j==2)swap2(a,i);
  if(j==4)swap4(a,i);
}

void readcn_(f1,a,ix,iy,n)      // array(i,j) form 0, as C
char f1[];
float a[];
int *ix,*iy,*n;
{
  FILE *fp;
  int i,j;
  int ip,jp,ii,jj=0;
  char b[36][80],f2[60];
  float c[2048],x;  
  for(i=0;i<60;i++){
    if(f1[i]<=' '){ f2[i]=0; break; }
    f2[i]=f1[i];
  }
  if((fp=fopen(f2,"r+b"))==NULL){ printf("open error\n"); return; }
  ip=0;
l10:
  ip++; if(ip>5){ fclose(fp); return; } 
  fread(b,36,80,fp);
  if(indexpos(b,"END     ",36)==36)goto l10;
  ip*=720;
  for(j=*iy;j<*iy+(*n);j++){
    if(j<0 || j>=2047)continue;
    jp=j*2048+(*ix)+ip;    
    fseek(fp,jp*4,0);
    fread(c,4,*n,fp);
    swap4(c,4*(*n));
    for(i=0;i<*n;i++){
      x=c[i];
      if(x<0.)x=0;
      a[jj++]=x;
    }
  }
  fclose(fp);
}

void readcc_(f1,a,ix,iy)
char f1[];
float a[101][101];
int *ix,*iy;
{
  FILE *fp;
  int i,j;
  int ip,jp,jx1,jx2,jy1,jy2,ii,jj;
  char b[36][80],f2[60];
  float c[101],x;  
  for(i=0;i<60;i++){
    if(f1[i]<=' '){ f2[i]=0; break; }
    f2[i]=f1[i];
  }
  if((fp=fopen(f2,"r+b"))==NULL){ printf("open error\n"); return; }
  ip=0;
l10:
  ip++; if(ip>5){ fclose(fp); return; } 
  fread(b,36,80,fp);
  if(indexpos(b,"END     ",36)==36)goto l10;
  ip*=720;
  for(i=0;i<101;i++)for(j=0;j<101;j++)a[i][j]=-1.;
  jx1=*ix-50; jx2=*ix+50;
  jy1=*iy-50; jy2=*iy+50;
  jj=101;
  for(j=jy1;j<=jy2;j++){
    jj--;
    if(j<0 || j>=2047)continue;
    jp=(j-1)*2048+jx1+ip-1;    
    fseek(fp,jp*4,0);
    fread(c,4,101,fp);
    swap4(c,404);
    for(i=0;i<101;i++){
      x=c[i];
      if(x<0.)x=0;
      a[jj][100-i]=x;
    }
  }
  fclose(fp);
}

void readc_(f1,a,ix,iy)
char f1[];
short a[101][101];
int *ix,*iy;
{
  FILE *fp;
  int i,j;
  int ip,jp,jx1,jx2,jy1,jy2,ii,jj;
  char b[36][80],f2[60];
  float c[101],x;  
  for(i=0;i<60;i++){
    if(f1[i]<=' '){ f2[i]=0; break; }
    f2[i]=f1[i];
  }
  if((fp=fopen(f2,"r+b"))==NULL){ printf("open error\n"); return; }
  ip=0;
l10:
  ip++; if(ip>5){ fclose(fp); return; } 
  fread(b,36,80,fp);
  if(indexpos(b,"END     ",36)==36)goto l10;
  ip*=720;
  for(i=0;i<101;i++)for(j=0;j<101;j++)a[i][j]=32767;
  jx1=*ix-50; jx2=*ix+50;
  jy1=*iy-50; jy2=*iy+50;
  jj=-1;
  for(j=jy1;j<=jy2;j++){
    jj++;
    if(j<0 || j>=2047)continue;
    jp=(j-1)*2048+jx1+ip-1;    
    fseek(fp,jp*4,0);
    fread(c,4,101,fp);
    swap4(c,404);
    for(i=0;i<101;i++){
      x=c[i]*0.25;
      if(x<0.)x=0;
      if(x>32766.)x=32766.;
      a[jj][i]=x;
    }
  }
  fclose(fp);
}

void readd_(f1,a,ix,iy,n1,n2)
char f1[];
short a[101][101];
int *ix,*iy,*n1,*n2;
{
  FILE *fp;
  int i,j;
  int ip,jp,jx1,jx2,jy1,jy2,ii,jj;
  char f2[40];
  short c[101],x;  
  for(i=0;i<60;i++){
    if(f1[i]<=' '){ f2[i]=0; break; }
    f2[i]=f1[i];
  }
  if((fp=fopen(f2,"r+b"))==NULL){ printf("open error\n"); return; }
  ip=3*1440;
  for(i=0;i<101;i++)for(j=0;j<101;j++)a[i][j]=32767;
  jx1=*ix-50; jx2=*ix+50;
  jy1=*iy-50; jy2=*iy+50;
  jj=101;
  for(j=jy1;j<=jy2;j++){
    jj--;
    if(j<0 || j>=(*n2)-1)continue;
    jp=(j-1)*(*n1)+jx1+ip-1;    
    fseek(fp,jp*2,0);
    fread(c,2,101,fp);
    swap2(c,202);
    for(i=0;i<101;i++){
      x=c[i];
      if(x<0)x=0;
      if(x>32766)x=32766;
      a[jj][i]=x;
    }
  }
  fclose(fp);
}

void wtn_(float a[], int *n1, int *n2, int *isign)
{
  int i1,i2,i3,k,n,nnew,nprev=1,nt,idim;
  float wksp[2048];
  n=*n1;  
  for (idim=0;idim<2;idim++) {
    nnew=n*nprev;
    for (i2=0;i2<(*n1)*(*n2);i2+=nnew) {
      for (i1=0;i1<nprev;i1++) {
        for (i3=i1+i2,k=0;k<n;k++,i3+=nprev) wksp[k]=a[i3];
        if (*isign >= 0) {
          for(nt=n;nt>=4;nt >>= 1) daub4(wksp,nt,*isign);
        } else {
          for(nt=4;nt<=n;nt <<= 1) daub4(wksp,nt,*isign);
        }
        for (i3=i1+i2,k=0;k<n;k++,i3+=nprev) a[i3]=wksp[k];
      }
    }
    nprev=nnew;
    n=*n2;
  }
}

daub4(float a[], int n, int isign)
{
  int nh,nh1,nh0,i,j;
  float wksp1[2048],C0,C1,C2,C3;
  C0=0.4829629131445341;
  C1=0.8365163037378079;
  C2=0.2241438680420134;
  C3=-.1294095225512604;
  nh=n>>1; nh1=nh+1; nh0=nh-1;
  if (isign >= 0) {
    for (i=0,j=0;j<n-3;j+=2,i++) {
    wksp1[i]    = C0*a[j]+C1*a[j+1]+C2*a[j+2]+C3*a[j+3];
    wksp1[i+nh] = C3*a[j]-C2*a[j+1]+C1*a[j+2]-C0*a[j+3];
  }
  wksp1[i]    = C0*a[n-2]+C1*a[n-1]+C2*a[0]+C3*a[1];
  wksp1[i+nh] = C3*a[n-2]-C2*a[n-1]+C1*a[0]-C0*a[1];
  } else {
    wksp1[0] = C2*a[nh0]+C1*a[n-1]+C0*a[0]+C3*a[nh];
    wksp1[1] = C3*a[nh0]-C0*a[n-1]+C1*a[0]-C2*a[nh];
    for (i=0,j=2; i<nh0; i++) {
      wksp1[j++] = C2*a[i]+C1*a[i+nh]+C0*a[i+1]+C3*a[i+nh1];
      wksp1[j++] = C3*a[i]-C0*a[i+nh]+C1*a[i+1]-C2*a[i+nh1];
    }
  }
  for (i=0;i<n;i++) a[i]=wksp1[i];
}