/****************************************/
/*	catio		valdes 8 30 82	*/
/*					*/
/*	Catalog io subroutines		*/
/****************************************/
#include	"focas1.h"


/*Open catalog file and return fd*/
FILE *
catopen (file, mode)
char   *file;
int     mode;
/*If mode = 0  Read*/
/*If mode = 1  Write*/
/*if mode = 2,3  Read / Write*/
{
    FILE *fd;

    switch (mode) {
    case 0: 
	if (strcmp (file, "STDIN") == 0)
	    fd = getstdin;
	else if ((fd = fopen (file, FOPEN_RO)) == NULL) {
	    focaserr (1, "catopen: Can't read catalog ", file);
	    return (NULL);
	}
	rdcathdr (fd, 0);
	return (fd);
    case 1: 
	if (strcmp (file, "STDOUT") == 0)
	    fd = stdout;
	else if ((fd = fopen (file, FOPEN_WO)) == NULL) {
	    c_delete (file);
	    if ((fd = fopen (file, FOPEN_WO)) == NULL) {
		focaserr (1, "catopen: Can't write to catalog ", file);
		return (NULL);
	    }
	}
	wtcathdr (fd, 0);
	return (fd);
    case 2: 
    case 3: 
	if (access(file,0)) {
	    if ((fd = fopen (file, FOPEN_WO)) == NULL) {
		focaserr (1, "catopen: Can't create catalog ", file);
		return (NULL);
	    }
	    wtcathdr (fd, 0);
	    fclose (fd);
	}
	if ((fd = fopen (file, FOPEN_RW)) == NULL) {
	    focaserr (1, "catopen: Can't read/write catalog ", file);
	    return (NULL);
	}
	rdcathdr (fd, 0);
	return (fd);
    default:
	focaserr (1, "catopen: Unknown mode", "");
	return (NULL);
    }
}

/*Read catalog header*/
rdcathdr (fd, pos)
FILE    *fd;
int     pos;
{
    int i, sz;
    char type[2];

    fseek (fd, 0L, pos);
    if (fread (type, 1, 2, fd) != 2)
	goto CATERR;
    cattype = type[0];

    switch (cattype) {
    case 1:
    case 2:
	fseek (fd, -1L, 1);
        sz = sizeof (struct sgparm);
        if (fread (&sp, sz, 1, fd) != 1)
	    goto CATERR;
        if (fread (&ncmmnts, 2, 1, fd) != 2)
	    goto CATERR;
        for (i = 0; i < ncmmnts; i++)
	    if (fread (comments[i], COLEN, 1, fd) != 1)
	        goto CATERR;
        sz = sizeof (struct flmcv);
        for (i = 0; i < sp.nden; i++)
        if (fread (&flmcv[i], sz, 1, fd) != 1)
	    goto CATERR;
        sz = sizeof (float) * sp.fx;
        for (i = 0; i < sp.fy; i++)
	    if (fread (dfltr[i], sz, 1, fd) != 1)
	        goto CATERR;
        sz = sizeof (float) * sp.nx;
        for (i = 0; i < sp.ny; i++)
	    if (fread (template[i], sz, 1, fd) != 1)
	        goto CATERR;
        sz = sizeof (struct crule);
        for (i = 0; i < sp.nrules; i++)
        if (fread (&crules[i], sz, 1, fd) != 1)
	    goto CATERR;
	break;
    case 3:
	fseek (fd, -1L, 1);
	if (mirdcathdr (fd) == 1)
	    goto CATERR;
	break;
    case 4:
	if (mirdcathdr (fd) == 1)
	    goto CATERR;
	break;
    default:
	focaserr (cattype, "rdcathdr: Unknown catalog type", "");
	return (1);
    }
    return (0);


CATERR:
    focaserr (5, "rdcathdr: Error reading catalog header", "");
    return (1);
}

/*Write catalog header */
wtcathdr (fd, pos)
FILE    *fd;
int     pos;
{
    char type[2];
    int i, sz;

    fseek (fd, 0L, pos);
    switch (cattype) {
    case 1:
    case 2:
	if (fwrite (&cattype, 1, 1, fd) != 1)
	    goto CATERR;
	sz = sizeof (struct sgparm);
        fwrite (&sp, sz, 1, fd);
        fwrite (&ncmmnts, 2, 1, fd);
        for (i = 0; i < ncmmnts; i++)
	    fwrite (comments[i], COLEN, 1, fd);
	sz = sizeof (struct flmcv);
        for (i = 0; i < sp.nden; i++)
            fwrite (&flmcv[i], sz, 1, fd);
	sz = sizeof (float) * sp.fx;
        for (i = 0; i < sp.fy; i++)
	    fwrite (dfltr[i], sz, 1, fd);
	sz = sizeof (float) * sp.nx;
        for (i = 0; i < sp.ny; i++)
	    fwrite (template[i], sz, 1, fd);
	sz = sizeof (struct crule);
        for (i = 0; i < sp.nrules; i++)
            fwrite (&crules[i], sz, 1, fd);
	break;
    case 3:
	if (fwrite (&cattype, 1, 1, fd) != 1)
	    goto CATERR;
	if (miwtcathdr (fd))
	    goto CATERR;
	break;
    case 4:
	type[0] = type[1] = 4;
	if (fwrite (type, 1, 2, fd) != 2)
	    goto CATERR;
	if (miwtcathdr (fd))
	    goto CATERR;
	break;
    }
    fseek (fd, 0L, 1);
    return;

CATERR:
    focaserr (6, "wtcathdr: Error in writing catalog header");
}

/*Size of catalog*/
long szcat (fd)
FILE    *fd;
{
    long    pos, end;

    pos = ftell (fd);
    fseek (fd, 0L, 2);
    end = ftell (fd);
    fseek (fd, pos, 0);
    return (end);
}

/*Size of header*/
long szcathdr () {
    long    size;

    switch (cattype) {
    case 1:
    case 2:
        size = sizeof (char);
        size += sizeof (struct sgparm);
        size += sizeof ncmmnts;
        size += ncmmnts * sizeof (char) * COLEN;
        size += sizeof (struct flmcv) * sp.nden;
        size += sizeof (float) * sp.fx * sp.fy;
        size += sizeof (float) * sp.nx * sp.ny;
        size += sp.nrules * sizeof (struct crule);
	return (size);
    case 3:
    case 4:
	return (miszcathdr());
    }
}

/*Number of objects in the catalog*/
long nobjs (fin)
FILE     *fin;
{
    long    szcat(), szcathdr(), minobjs();

    switch (cattype) {
    case 1: 
	return ((long) (szcat(fin) - szcathdr()) / (sizeof (struct objrec)));
    case 2: 
	return ((long) (szcat(fin) - szcathdr()) / (sizeof (struct shortrec)));
    case 3: 
    case 4: 
	return (minobjs (fin));
    }
}

long szobjrec ()
{
    switch (cattype) {
    case 1:
	return (sizeof (struct objrec));
    case 2:
	return (sizeof (struct shortrec));
    case 3:
    case 4:
	return (miszobjrec());
    }
}

/*Read an object record with possible offset*/
/*Return 1 at eof, -1 if seek before first object*/
rdcatob (fd, offset, ob)
FILE    *fd;
long    offset;
struct objrec  *ob;
{
    long    pos, szobjrec;
    struct shortrec sr;

    switch (cattype) {
    case 1: 
	pos = ftell (fd);
	szobjrec = sizeof (struct objrec);
	offset *= szobjrec;
	if (offset == 0) {
	    if (fread (ob, szobjrec, 1, fd) < 1)
		return (1);
	} else if (offset > 0) {
	    if (pos + offset >= szcat (fd))
		return (1);
	    fseek (fd, offset, 1);
	    if (fread (ob, szobjrec, 1, fd) < 1)
		return (1);
	} else {
	    if (pos + offset < szcathdr ())
		return (-1);
	    fseek (fd, offset, 1);
	    if (fread (ob, szobjrec, 1, fd) < 1)
		return (1);
	}
	break;
    case 2: 
	pos = ftell (fd);
	szobjrec = sizeof (struct shortrec);
	offset *= sizeof (struct shortrec);
	if (offset == 0) {
	    if (fread (&sr, szobjrec, 1, fd) < 1)
		return (1);
	} else if (offset > 0) {
	    if (pos + offset >= szcat (fd))
		return (1);
	    fseek (fd, offset, 1);
	    if (fread (&sr, szobjrec, 1, fd) < 1)
		return (1);
	} else {
	    if (pos + offset < szcathdr ())
		return (-1);
	    fseek (fd, offset, 1);
	    if (fread (&sr, szobjrec, 1, fd) < 1)
		return (1);
	}
	srtolr (&sr, ob);
	break;
    case 3: 
    case 4: 
	if (mirdcatob (fd, offset, ob) == 1)
	    return (1);
	break;
    }
    return (0);
}

wtcat (ob, fd)
struct objrec  *ob;
FILE    *fd;
{
    struct shortrec sr;

    fseek (fd, 0L, 1);
    switch (cattype) {
    case 1: 
	if (fwrite (ob, sizeof (struct objrec), 1, fd) != 1) {
	    focaserr (1, "wtcat: Error on write", "");
	    return;
	}
	break;
    case 2: 
	lrtosr (ob, &sr);
	if (fwrite (&sr, sizeof (struct shortrec), 1, fd) != 1) {
	    focaserr (1, "wtcat: Error on write", "");
	    return;
	}
	break;
    case 3: 
    case 4: 
	if (miwtcat (fd, ob) == 1) {
	    focaserr (1, "wtcat: Error on write", "");
	    return;
	}
	break;
    }
    fseek (fd, 0L, 1);
}

/*search catalog for particular entry
rdentry(entnum, subent, ob, fd)
  short	entnum;
  long	subent;
  struct objrec	*ob;
  FILE	*fd;
  {
  int	i;

  fprintf(stderr, "rdentry: search for %d.%-d\n", entnum, subent);
  fseek(fd, 0L, 0);
  rdcathdr(fd, 0);
  for (;;) {
    rdcatob(fd, 0, ob);
    if (feof(fd)) {
      fprintf(stderr, "rdentry: search failed\n");
      return(1);
      }
    if ((ob->entnum == entnum) && (ob->subent == subent)) return(0);
    }
  }
*/

/*Return 1 if entry 1 is parent, 2 if entry 2 is parent, and 0 if unrelated*/
subentry (e1, s1, e2, s2)
short   e1, e2;			/* entry numbers */
long    s1, s2;			/* subentry numbers */
{
    int     i;

    if (e1 != e2)
	return (0);
    if ((s1 == 0) && (s2 != 0))
	return (1);
    if ((s2 == 0) && (s1 != 0))
	return (2);
    i = s2;
    for (i = s2; i > 0; i /= 10)
	if (i == s1)
	    return (1);
    for (i = s1; i > 0; i /= 10)
	if (i == s2)
	    return (2);
    return (0);
}

/*Convert short record to long record*/
srtolr (sr, lr)
struct objrec  *lr;
struct shortrec *sr;
{

    lr -> entnum = sr -> entnum;
    lr -> subent = sr -> subent;
    strcpy (lr -> class, sr -> class);
    lr -> eflgs = sr -> eflgs;
    lr -> xc = sr -> xc;
    lr -> yc = sr -> yc;
    lr -> mag = sp.magfst - sr -> magi / 100.;
    if (sr -> magi > -10000)
	lr -> Li = pow (10., (sr -> magi / 250.));
    else
	lr -> Li = 0.;
    if (sr -> magc > -10000)
	lr -> Lc = pow (10., (sr -> magc / 250.));
    else
	lr -> Lc = 0.;
    if (sr -> magfca > -10000)
	lr -> Lfca = pow (10., (sr -> magfca / 250.));
    else
	lr -> Lfca = 0.;
    if (sr -> magtotal > -10000)
	lr -> Ltotal = pow (10., (sr -> magtotal / 250.));
    else
	lr -> Ltotal = 0.;
    lr -> icx = sr -> icx;
    lr -> icy = sr -> icy;
    lr -> ixx = sr -> ixx;
    lr -> iyy = sr -> iyy;
    lr -> ixy = sr -> ixy;
    lr -> ir1 = sr -> ir1;
    lr -> cx = sr -> cx;
    lr -> cy = sr -> cy;
    lr -> xx = sr -> xx;
    lr -> yy = sr -> yy;
    lr -> xy = sr -> xy;
    lr -> r1 = sr -> r1;
    lr -> scale = sr -> scale / 100.;
    lr -> frac = sr -> frac / 100.;
    lr -> area = sr -> area;	/* JFJarvis 3/11/85 */
}

/*Convert long record to short record*/
lrtosr (lr, sr)
struct objrec  *lr;
struct shortrec *sr;
{
    float   intmag ();

    sr -> entnum = lr -> entnum;
    sr -> subent = lr -> subent;
    strcpy (sr -> class, lr -> class);
    sr -> eflgs = lr -> eflgs;
    sr -> xc = lr -> xc;
    sr -> yc = lr -> yc;
    if (lr -> Li > 0.)
	sr -> magi = 250 * log10 (lr -> Li);
    else
	sr -> magi = -10000;
    if (lr -> Lc > 0.)
	sr -> magc = 250 * log10 (lr -> Lc);
    else
	sr -> magc = -10000;
    if (lr -> Lfca > 0.)
	sr -> magfca = 250 * log10 (lr -> Lfca);
    else
	sr -> magfca = -10000;
    if (lr -> Ltotal > 0.)
	sr -> magtotal = 250 * log10 (lr -> Ltotal);
    else
	sr -> magtotal = -10000;
    sr -> icx = lr -> icx;
    sr -> icy = lr -> icy;
    sr -> ixx = lr -> ixx;
    sr -> iyy = lr -> iyy;
    sr -> ixy = lr -> ixy;
    sr -> ir1 = lr -> ir1;
    sr -> cx = lr -> cx;
    sr -> cy = lr -> cy;
    sr -> xx = lr -> xx;
    sr -> yy = lr -> yy;
    sr -> xy = lr -> xy;
    sr -> r1 = lr -> r1;
    sr -> scale = 100 * lr -> scale;
    sr -> frac = 100 * lr -> frac;
    sr -> area = lr -> area;	/* JFJarvis 3/11/85 */
}

/*Print transformation matrix*/
prxfrm () {
    int     i, j;

    printf ("\nTransformation matrix\n\n");
    for (i = 0; i < 3; i++) {
	for (j = 0; j < 3; j++)
	    printf ("%12g ", sp.xfrm[i][j]);
	printf ("\n");
    }
}

/*Print detection filter*/
prdflt () {
    int     k, l;

    printf ("\nDetection filter\n\n");
    if (sp.pflags & FLTR)
	printf ("  Builtin filter\n");
    for (k = 0; k < sp.fy; k++) {
	for (l = 0; l < sp.fx; l++)
	    printf ("%5g ", dfltr[k][l]);
	printf ("\n");
    }
}

/*Print class rules*/
prclssr () {
    int     k;

    printf ("\nClassification rules:\n\n");
    for (k = 0; k < sp.nrules; k++)
	printf ("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %5s\n",
		crules[k].mag1, crules[k].mag2, crules[k].scl1, crules[k].scl2,
		crules[k].frc1, crules[k].frc2, crules[k].class);
}

/*Print psf*/
prpsf () {
    int     k, l;
    float   norm;

    printf ("\nPoint spread function template\n\n");
    norm = template[sp.ny / 2][sp.nx / 2];
    for (k = 0; k < sp.ny; k++) {
	for (l = 0; l < sp.nx; l++)
	    printf ("%3.0f ", 100.* template[k][l] / norm);
	printf ("\n");
    }
}

/*Print header*/
prhdr () {
    int     k;
    float   corners[4][2], recorn[4][2];

    printf ("Field image file: %s\n", sp.ptfl);
    printf ("Field name: %s\n", sp.pltnm);
    printf ("Field epoch: %s\n", sp.pltepoch);
    printf ("Field passband: %s\n", sp.band);
    printf ("Field coordinates: %.6g %.6g\n", sp.pltra, sp.pltdec);
    printf ("Image exposure or integration: %d\n", sp.pltexp);
    printf ("Observer: %s\n", sp.observer);
    printf ("Origin: %s\n", sp.origin);
    printf ("Field size: %d %d\n", sp.kpts, sp.kscan);
    printf ("Field scan origin: %d %d\n",
	(int) lpcoords (0., 1, 1), (int) lpcoords (0., 2, 1));
    corners[0][0] = lpcoords (0., 1, 1);
    corners[0][1] = lpcoords (0., 2, 1);
    corners[1][0] = lpcoords (sp.kpts-1., 1, 1);
    corners[1][1] = lpcoords (0., 2, 1);
    corners[2][0] = lpcoords (0., 1, 1);
    corners[2][1] = lpcoords (sp.kscan-1., 2, 1);
    corners[3][0] = lpcoords (sp.kpts-1., 1, 1);
    corners[3][1] = lpcoords (sp.kscan-1., 2, 1);
    printf ("Field corners in reference coordinates:\n");
    for (k = 0; k < 4; k++) {
	/*
	coords (corners[k], recorn[k]);
	printf (" (%.6g, %.6g)", recorn[k][0], recorn[k][1]);
	*/
	printf (" (%.6g, %.6g)", corners[k][0], corners[k][1]);
    }
    printf ("\n");
    printf ("Magnitude zero point: %f\n", sp.magfst);
    printf ("Catalog magnitude limit: %4.2f\n", sp.maglim);
    printf ("Radius of fixed circular aperature: %f\n", sp.rfca);
    printf ("Sigma of sky: %.6g\n", sp.skyhw);
    printf ("Sigma above sky for detection: %f\n", sp.thrln0);
    printf ("Sigma below sky for detection: %f\n", sp.thrln1);
    printf ("Minimum area for detection: %d\n", sp.minarea);
    printf ("Significance for evaluation and splitting: %.6g\n", sp.sig);
    printf ("Area file: %12s\n", sp.arfl);
    printf ("\nComments:\n\n");
    for (k = 0; k < sp.comment % ncmmnts; k++)
	printf ("%s\n", comments[k]);
}

/*Print intensity relation*/
printen () {
    int     k;

    printf ("\nInstrumental intensity relation:\n");
    if (sp.pflags & LNR) {
	printf ("\n  Linear\n");
	return;
    }
    else
	if (sp.pflags & DLUT)
	    printf ("\n  Direct Look Up Table\n");
	else
	    if (sp.pflags & INTP)
		printf ("\n  Interpolation Look Up Table\n");
	    else
		printf ("\n  Function evaluation\n");
/*    if ((sp.pflags & DLUT) || (sp.pflags & INTP)) {
	for (k = 0; k < sp.nden; k++)
	    printf ("D-TO-I:  %3d  %.6g\n", flmcv[k].d, flmcv[k].i);
    } */
}

/* Convert instrumental luminosity to magnitude */
float
        intmag (I)
float   I;
{
    if (I > 0)
	return (sp.magfst - 2.5 * log10 (I));
    else
	if (I == 0.)
	    return (100.);
	else
	    return (sp.magfst - 2.5 * log10 (-I));
}