diff options
Diffstat (limited to 'noao/nproto/ace/catio.x')
-rw-r--r-- | noao/nproto/ace/catio.x | 931 |
1 files changed, 931 insertions, 0 deletions
diff --git a/noao/nproto/ace/catio.x b/noao/nproto/ace/catio.x new file mode 100644 index 00000000..1fbae947 --- /dev/null +++ b/noao/nproto/ace/catio.x @@ -0,0 +1,931 @@ +include <imset.h> +#include <tbset.h> +define TBL_NROWS 0 +include <math.h> +include "ace.h" +include "cat.h" +include "objs.h" + + +# CATOPEN -- Open a catalog. +# This may be used just to allocate the structure or to actually open +# a catalog file. It does not read the objects. Use catrobjs. + +procedure catopen (cat, input, output, catdef) + +pointer cat #U Catalog structure +char input[ARB] #I Input catalog name +char output[ARB] #I Output catalog name +char catdef[ARB] #I Catalog definition file + +pointer tbl + +bool streq() +pointer tbtopn() + +begin + if (cat == NULL) + call calloc (cat, CAT_LEN, TY_STRUCT) + + if (input[1] == EOS && output[1] == EOS) + return + + if (streq (input, output)) { # READ_WRITE + call calloc (tbl, TBL_LEN, TY_STRUCT) + CAT_INTBL(cat) = tbl + CAT_OUTTBL(cat) = tbl + + TBL_TP(tbl) = tbtopn (input, READ_WRITE, 0) + call catdefine (tbl, READ_ONLY, catdef) + call catrhdr (cat) + } else if (output[1] == EOS) { # READ_ONLY + call calloc (tbl, TBL_LEN, TY_STRUCT) + CAT_INTBL(cat) = tbl + CAT_OUTTBL(cat) = NULL + + TBL_TP(tbl) = tbtopn (input, READ_ONLY, 0) + call catdefine (tbl, READ_ONLY, catdef) + call catrhdr (cat) + } else if (input[1] == EOS) { # NEW_FILE + call calloc (tbl, TBL_LEN, TY_STRUCT) + CAT_INTBL(cat) = NULL + CAT_OUTTBL(cat) = tbl + + TBL_TP(tbl) = tbtopn (output, NEW_FILE, 0) + call catdefine (tbl, NEW_FILE, catdef) + } else { # NEW_COPY + call calloc (tbl, TBL_LEN, TY_STRUCT) + CAT_INTBL(cat) = tbl + + TBL_TP(tbl) = tbtopn (input, READ_ONLY, 0) + call catdefine (tbl, NEW_COPY, catdef) + call catrhdr (cat) + + call calloc (tbl, TBL_LEN, TY_STRUCT) + CAT_OUTTBL(cat) = tbl + TBL_TP(tbl) = tbtopn (output, NEW_COPY, TBL_TP(CAT_INTBL(cat))) + call catdefine (tbl, NEW_COPY, catdef) + } +end + + +procedure catcreate (cat) + +pointer cat #I Catalog structure + +pointer tbl, tp + +begin + if (cat == NULL) + return + tbl = CAT_OUTTBL(cat) + if (tbl == NULL) + return + tp = TBL_TP(tbl) + if (tp == NULL) + return + if (CAT_INTBL(cat) != NULL) { + if (tp == TBL_TP(CAT_INTBL(cat))) + return + } + call tbtcre (tp) +end + + +# CATCLOSE -- Close a catalog. + +procedure catclose (cat) + +pointer cat #I Catalog pointer + +int i +pointer tbl, objs + +begin + if (cat == NULL) + return + + tbl = CAT_INTBL(cat) + if (tbl != NULL) { + if (TBL_STP(tbl) != NULL) + call stclose (TBL_STP(tbl)) + if (tbl == CAT_OUTTBL(cat)) + CAT_OUTTBL(cat) = NULL + call tbtclo (TBL_TP(tbl)) + } + tbl = CAT_OUTTBL(cat) + if (tbl != NULL) { + if (TBL_STP(tbl) != NULL) + call stclose (TBL_STP(tbl)) + call tbtclo (TBL_TP(tbl)) + } + + objs = CAT_OBJS(cat) + if (objs != NULL) { + do i = 0, CAT_NUMMAX(cat)-1 + call mfree (Memi[objs+i], TY_STRUCT) + } + + call mfree (CAT_APFLUX(cat), TY_REAL) + call mfree (CAT_OBJS(cat), TY_POINTER) + call mfree (CAT_INTBL(cat), TY_STRUCT) + call mfree (CAT_OUTTBL(cat), TY_STRUCT) + call mfree (CAT_HDR(cat), TY_STRUCT) + call mfree (cat, TY_STRUCT) +end + + +# CATGETS -- Get a string parameter from the catalog header. + +procedure catgets (cat, param, value, maxchar) + +pointer cat #I Catalog pointer +char param[ARB] #I Parameter to get +char value[ARB] #O Returned value +int maxchar #I Maximum characters in value + +int i, strdic() + +begin + value[1] = EOS + + if (cat == NULL) + return + + i = strdic (param, CAT_STR(cat), CAT_SZSTR, CATPARAMS) + switch (i) { + case 1: + if (CAT_HDR(cat) == NULL) + i = 0 + else + call strcpy (HDR_IMAGE(CAT_HDR(cat)), value, maxchar) + case 2: + if (CAT_HDR(cat) == NULL) + i = 0 + else + call strcpy (HDR_MASK(CAT_HDR(cat)), value, maxchar) + case 3: + call strcpy (CAT_OBJID(cat), value, maxchar) + case 4: + call strcpy (CAT_CATALOG(cat), value, maxchar) + default: + call sprintf (CAT_STR(cat), CAT_SZSTR, + "catgets: unknown catalog parameter `%s'") + call pargstr (param) + call error (1, CAT_STR(cat)) + } + + if (i == 0) { + call sprintf (CAT_STR(cat), CAT_SZSTR, + "catgets: parameter `%s' not found") + call pargstr (param) + call error (1, CAT_STR(cat)) + } +end + + +procedure catgeti (cat, param, value) + +pointer cat #I Catalog pointer +char param[ARB] #I Parameter to get +int value #O Returned value + +int i, strdic() + +begin + value = INDEFI + + if (cat == NULL) + return + + i = strdic (param, CAT_STR(cat), CAT_SZSTR, CATPARAMS) + switch (i) { + case 5: + value = CAT_NOBJS(cat) + default: + call sprintf (CAT_STR(cat), CAT_SZSTR, + "catgeti: unknown catalog parameter `%s'") + call pargstr (param) + call error (1, CAT_STR(cat)) + } +end + + +procedure catputs (cat, param, value) + +pointer cat #I Catalog pointer +char param[ARB] #I Parameter to get +char value[ARB] #I Value + +int i, strdic() + +begin + if (cat == NULL) + return + + i = strdic (param, CAT_STR(cat), CAT_SZSTR, CATPARAMS) + switch (i) { + case 0: + call sprintf (CAT_STR(cat), CAT_SZSTR, + "catgets: unknown catalog parameter `%s'") + call pargstr (param) + call error (1, CAT_STR(cat)) + case 1: + if (CAT_HDR(cat) == NULL) + call calloc (CAT_HDR(cat), HDR_LEN, TY_STRUCT) + call strcpy (value, HDR_IMAGE(CAT_HDR(cat)), HDR_SZFNAME) + case 2: + if (CAT_HDR(cat) == NULL) + call calloc (CAT_HDR(cat), HDR_LEN, TY_STRUCT) + call strcpy (value, HDR_MASK(CAT_HDR(cat)), HDR_SZFNAME) + case 3: + call strcpy (value, CAT_OBJID(cat), CAT_SZSTR) + case 4: + call strcpy (value, CAT_CATALOG(cat), CAT_SZSTR) + } +end + + +procedure catputr (cat, param, value) + +pointer cat #I Catalog pointer +char param[ARB] #I Parameter to get +real value #I Value + +int i, strdic() + +begin + if (cat == NULL) + return + + i = strdic (param, CAT_STR(cat), CAT_SZSTR, CATPARAMS) + switch (i) { + case 6: + if (CAT_HDR(cat) == NULL) + call calloc (CAT_HDR(cat), HDR_LEN, TY_STRUCT) + HDR_MAGZERO(CAT_HDR(cat)) = value + default: + call sprintf (CAT_STR(cat), CAT_SZSTR, + "catgetr: unknown catalog parameter `%s'") + call pargstr (param) + call error (1, CAT_STR(cat)) + } +end + + +procedure catrobjs (cat, filt) + +pointer cat #I Catalog pointer +char filt[ARB] #I Filter string + +int i, num, nrows, nobjs, nummax, nalloc, tbpsta() +pointer tbl, tp, objs, obj +bool filter() + +begin + if (cat == NULL) + return + + tbl = CAT_INTBL(cat) + if (tbl == NULL) + return + tp = TBL_TP(tbl) + + nrows = tbpsta (tp, TBL_NROWS) + nalloc = nrows + NUMSTART - 1 + call calloc (objs, nalloc, TY_POINTER) + + nobjs = 0 + nummax = 0 + obj = NULL + do i = 1, nrows { + call catrobj (cat, obj, i) + if (!filter (obj, filt)) + next + num = OBJ_NUM(obj) + if (num > nalloc) { + nalloc = nalloc + 1000 + call realloc (objs, nalloc, TY_POINTER) + call aclri (Memi[objs+nalloc-1000], 1000) + } + if (Memi[objs+num-1] == NULL) + nobjs = nobjs + 1 + nummax = max (num, nummax) + Memi[objs+num-1] = obj + obj = NULL + } + + CAT_OBJS(cat) = objs + CAT_NOBJS(cat) = nobjs + CAT_NUMMAX(cat) = nummax +end + + +procedure catrobj (cat, obj, row) + +pointer cat #I Catalog pointer +pointer obj #U Object pointer +int row #I Table row + +int id, type, ori() +pointer tbl, tp, stp, sym, cdef, sthead(), stnext() + +begin + if (cat == NULL) + return + + tbl = CAT_INTBL(cat) + if (tbl == NULL) + return + + tp = TBL_TP(tbl) + stp = TBL_STP(tbl) + + if (obj == NULL) + call calloc (obj, OBJ_LEN, TY_STRUCT) + + for (sym=sthead(stp); sym!=NULL; sym=stnext(stp,sym)) { + id = ENTRY_ID(sym) + if (id > 1000 || id == ID_APFLUX) + next + switch (id) { + case ID_FLAGS: + OBJ_FLAGS(obj) = 0 + ifnoerr (call tbegtt (tp, cdef, row, CAT_STR(cat), CAT_SZSTR)) { + if (Memc[CAT_STRPTR(cat)] == 'M') + SETFLAG(obj,OBJ_SPLIT) + } + next + } + + type = ENTRY_TYPE(sym) + cdef = ENTRY_CDEF(sym) + switch (type) { + case TY_INT: + iferr (call tbegti (tp, cdef, row, OBJI(obj,id))) + OBJI(obj,id) = INDEFI + case TY_REAL: + iferr (call tbegtr (tp, cdef, row, OBJR(obj,id))) + OBJR(obj,id) = INDEFR + case TY_DOUBLE: + iferr (call tbegtd (tp, cdef, row, OBJD(obj,id))) + OBJD(obj,id) = INDEFD + default: + iferr (call tbegtt (tp, cdef, row, OBJC(obj,id), -type)) + OBJC(obj,id) = EOS + } + } + + OBJ_ROW(obj) = row +end + + +procedure catwobj (cat, obj, row) + +pointer cat #I Catalog pointer +pointer obj #I Object pointer +int row #I Table row + +int ival +real rval +double dval +pointer sval + +int id, type, func, napr, andi() +real magzero, a, b, theta, elong, ellip, r, cxx, cyy, cxy +real aerr, berr, thetaerr, cxxerr, cyyerr, cxyerr +bool doshape +pointer tbl, tp, stp, sym, cdef, sthead(), stnext() + +begin + if (obj == NULL) + return + + tbl = CAT_OUTTBL(cat) + if (tbl == NULL) + return + tp = TBL_TP(tbl) + stp = TBL_STP(tbl) + + #call sprintf (CAT_STR(cat), CAT_SZSTR, "%s-%d") + # if (OBJ_OBJID(obj) != NULL) + # call pargstr (Memc[OBJ_OBJID(obj)]) + # else + # call pargstr (CAT_OBJID(cat)) + # call pargi (OBJ_NUM(obj)) + #call tbeptt (tp, TBL_BJID(tbl), row, CAT_STR(cat)) + #call tbeptt (tp, TBL_CLASS(tbl), row, OBJ_CLASS(obj)) + + magzero = CAT_MAGZERO(cat) + if (IS_INDEFR(magzero)) + magzero = 0. + sval = CAT_STRPTR(cat) + napr = 0 + doshape = false + for (sym=sthead(stp); sym!=NULL; sym=stnext(stp,sym)) { + id = ENTRY_ID(sym) + func = ENTRY_FUNC(sym) + type = ENTRY_TYPE(sym) + cdef = ENTRY_CDEF(sym) + if (id > 1000) { + switch (id) { + case ID_A, ID_B, ID_THETA, ID_ELONG, ID_ELLIP, ID_R, ID_CXX, + ID_CYY, ID_CXY: + if (!doshape) { + call catshape (obj, a, b, theta, elong, ellip, r, + cxx, cyy, cxy, aerr, berr, thetaerr, cxxerr, + cyyerr, cxyerr) + doshape = true + } + switch (id) { + case ID_A: + rval = a + case ID_B: + rval = b + case ID_THETA: + rval = theta + case ID_ELONG: + rval = elong + case ID_ELLIP: + rval = ellip + case ID_R: + rval = r + case ID_CXX: + rval = cxx + case ID_CYY: + rval = cyy + case ID_CXY: + rval = cxy + } + case ID_FLUXERR, ID_XERR, ID_YERR: + switch (id) { + case ID_FLUXERR: + rval = OBJ_FLUXVAR(obj) + case ID_XERR: + rval = OBJ_XVAR(obj) + case ID_YERR: + rval = OBJ_YVAR(obj) + } + if (IS_INDEFR(rval) || rval < 0.) + rval = INDEFR + else + rval = sqrt (rval) + case ID_AERR, ID_BERR, ID_THETAERR, ID_CXXERR, ID_CYYERR, + ID_CXYERR: + if (!doshape) { + call catshape (obj, a, b, theta, elong, ellip, r, + cxx, cyy, cxy, aerr, berr, thetaerr, cxxerr, + cyyerr, cxyerr) + doshape = true + } + switch (id) { + case ID_AERR: + rval = aerr + case ID_BERR: + rval = aerr + case ID_THETAERR: + rval = aerr + case ID_CXXERR: + rval = aerr + case ID_CYYERR: + rval = aerr + case ID_CXYERR: + rval = aerr + } + } + } else if (id == ID_FLAGS) { + if (SPLIT(obj)) + call strcpy ("M", Memc[sval], SZ_LINE) + else + call strcpy ("-", Memc[sval], SZ_LINE) + } else if (id == ID_APFLUX) { + if (OBJ_APFLUX(obj) == NULL) + rval = INDEFR + else { + rval = Memr[OBJ_APFLUX(obj)+napr] + napr = napr + 1 + } + } else { + switch (type) { + case TY_INT: + ival = OBJI(obj,id) + case TY_REAL: + rval = OBJR(obj,id) + case TY_DOUBLE: + dval = OBJD(obj,id) + default: + call strcpy (OBJC(obj,id), Memc[sval], SZ_LINE) + } + } + + # Apply function. + if (func > 0) { + if (ENTRY_CTYPE(sym) != type) { + # For now all function types are real. + switch (type) { + case TY_INT: + rval = ival + case TY_DOUBLE: + rval = dval + } + } + type = ENTRY_CTYPE(sym) + switch (func) { + case FUNC_MAG: + if (!IS_INDEFR(rval)) { + if (rval <= 0.) + rval = INDEFR + else + rval = -2.5 * log10 (rval) + magzero + } + } + } + + # Write to catalog. + switch (type) { + case TY_INT: + call tbepti (tp, cdef, row, ival) + case TY_REAL: + call tbeptr (tp, cdef, row, rval) + case TY_DOUBLE: + call tbeptd (tp, cdef, row, dval) + default: + call tbeptt (tp, cdef, row, Memc[sval]) + } + } + OBJ_ROW(obj) = row +end + + +# CATWCS -- Set catalog WCS information. + +procedure catwcs (cat, im) + +pointer cat #I Catalog pointer +pointer im #I IMIO pointer + +int i +pointer sp, axtype, label, units, format +pointer mw, tbl, tp, stp, sym, cdef + +bool streq() +pointer mw_openim(), sthead(), stnext(), stname() +errchk mw_openim + +begin + if (cat == NULL) + return + if (CAT_OUTTBL(cat) == NULL) + return + + call smark (sp) + call salloc (axtype, SZ_FNAME, TY_CHAR) + call salloc (label, SZ_FNAME, TY_CHAR) + call salloc (units, SZ_FNAME, TY_CHAR) + call salloc (format, SZ_FNAME, TY_CHAR) + + tbl = CAT_OUTTBL(cat) + tp = TBL_TP(tbl) + stp = TBL_STP(tbl) + + mw = mw_openim (im) + do i = 1, 2 { + iferr (call mw_gwattrs (mw, i, "axtype", Memc[axtype], SZ_FNAME)) + Memc[axtype] = EOS + iferr (call mw_gwattrs (mw, i, "label", Memc[label], SZ_FNAME)) { + if (streq (Memc[axtype], "ra")) + call strcpy ("RA", Memc[label], SZ_FNAME) + else if (streq (Memc[axtype], "dec")) + call strcpy ("DEC", Memc[label], SZ_FNAME) + else + Memc[label] = EOS + } + iferr (call mw_gwattrs (mw, i, "units", Memc[units], SZ_FNAME)) { + if (streq (Memc[axtype], "ra") || streq (Memc[axtype], "dec")) + call strcpy ("deg", Memc[units], SZ_FNAME) + else + Memc[units] = EOS + } + iferr (call mw_gwattrs (mw, i, "format", Memc[format], SZ_FNAME)) { + if (streq (Memc[axtype], "ra")) + call strcpy ("%.2H", Memc[format], SZ_FNAME) + else if (streq (Memc[axtype], "dec")) + call strcpy ("%.1h", Memc[format], SZ_FNAME) + else + Memc[format] = EOS + } + + if (i == 1) { + for (sym=sthead(stp); sym!=NULL; sym=stnext(stp,sym)) { + if (ENTRY_ID(sym) != ID_WX) + next + if (!(streq (Memc[stname(stp,sym)], "WX") || + streq (Memc[stname(stp,sym)], "wx"))) + Memc[label] = EOS + break + } + } else { + for (sym=sthead(stp); sym!=NULL; sym=stnext(stp,sym)) { + if (ENTRY_ID(sym) != ID_WY) + next + if (!(streq (Memc[stname(stp,sym)], "WY") || + streq (Memc[stname(stp,sym)], "wy"))) + Memc[label] = EOS + break + } + } + + if (sym != NULL) { + cdef = ENTRY_CDEF(sym) + if (Memc[label] != EOS) + call tbcnam (tp, cdef, Memc[label]) + if (Memc[units] != EOS) + call tbcnit (tp, cdef, Memc[units]) + if (Memc[format] != EOS) + call tbcfmt (tp, cdef, Memc[format]) + } + } + call mw_close (mw) + + call sfree (sp) +end + + +procedure catrhdr (cat) + +pointer cat #I Catalog pointer + +pointer tp, hdr + +begin + if (cat == NULL) + return + + if (CAT_HDR(cat) != NULL) + call mfree (CAT_HDR(cat), TY_STRUCT) + if (CAT_INTBL(cat) == NULL) + return + tp = TBL_TP(CAT_INTBL(cat)) + + call calloc (CAT_HDR(cat), HDR_LEN, TY_STRUCT) + hdr = CAT_HDR(cat) + + iferr (call tbhgtt (tp, "IMAGE", HDR_IMAGE(hdr), HDR_SZFNAME)) + HDR_IMAGE(hdr) = EOS + iferr (call tbhgtt (tp, "MASK", HDR_MASK(hdr), HDR_SZFNAME)) + HDR_MASK(hdr) = EOS + iferr (call tbhgtr (tp, "MAGZERO", HDR_MAGZERO(hdr))) + HDR_MAGZERO(hdr) = INDEFR +end + + +procedure catwhdr (cat, im) + +pointer cat #I Catalog pointer +pointer im #I Image pointer + +pointer tp, hdr + +begin + if (cat == NULL) + return + + tp = CAT_OUTTBL(cat) + hdr = CAT_HDR(cat) + if (tp == NULL || hdr == NULL) + return + tp = TBL_TP(tp) + + if (HDR_IMAGE(hdr) != EOS) + call tbhadt (tp, "IMAGE", HDR_IMAGE(hdr)) + if (HDR_MASK(hdr) != EOS) + call tbhadt (tp, "MASK", HDR_MASK(hdr)) + if (!IS_INDEFR(HDR_MAGZERO(hdr))) + call tbhadr (tp, "MAGZERO", HDR_MAGZERO(hdr)) +end + + +procedure catwobjs (cat) + +pointer cat #I Catalog pointer + +int i, j +pointer objs, obj + +begin + if (cat == NULL) + return + if (CAT_OUTTBL(cat) == NULL) + return + if (CAT_OBJS(cat) == NULL) + return + + objs = CAT_OBJS(cat) + j = 0 + do i = 1, CAT_NUMMAX(cat) { + obj = Memi[objs+i-1] + if (obj == NULL) + next + j = j + 1 + call catwobj (cat, obj, j) + } +end + + + +procedure catdump (cat) + +pointer cat #I Catalog pointer + +int i +pointer objs, obj + +begin + if (CAT_OBJS(cat) == NULL) + return + + objs = CAT_OBJS(cat) + do i = 1, CAT_NUMMAX(cat) { + obj = Memi[objs+i-1] + if (obj == NULL) + next + + call printf ("%d %d %g %g\n") + call pargi (OBJ_NUM(obj)) + call pargi (OBJ_NPIX(obj)) + call pargr (OBJ_XAP(obj)) + call pargr (OBJ_YAP(obj)) + } +end + + +# CATGOBJ -- Get object given the object number. +# +# Currently this relies on the object pointer array being indexed by +# object number. + +pointer procedure catgobj (cat, num) + +pointer cat #I Catalog +int num #I Object number + +begin + return (Memi[CAT_OBJS(cat)+num-1]) +end + + +# These currently work on the object number but eventually there will be +# an array of indices to allow traversing the objects in some sorted order. + +pointer procedure cathead (cat) + +pointer cat #I Catalog pointer + +int i +pointer objs, obj + +begin + objs = CAT_OBJS(cat) + do i = 0, CAT_NUMMAX(cat)-1 { + obj = Memi[objs+i] + if (obj != NULL) + return (obj) + } + return (NULL) +end + + +pointer procedure catnext (cat, obj) + +pointer cat #I Catalog pointer +pointer obj #I Input object pointer + +int i +pointer objs, objnext + +begin + if (obj == NULL) + return (NULL) + + objs = CAT_OBJS(cat) + do i = OBJ_NUM(obj), CAT_NUMMAX(cat)-1 { + objnext = Memi[objs+i] + if (objnext != NULL) + return (objnext) + } + return (NULL) +end + + +procedure catshape (obj, a, b, theta, elong, ellip, r, cxx, cyy, cxy, + aerr, berr, thetaerr, cxxerr, cyyerr, cxyerr) + +pointer obj #I Object structure +real a #O Semimajor axis based on second moments +real b #O Semiminor axis based on second moments +real theta #O Position angle based on second moments +real elong #O Elongation (A/B) +real ellip #O Ellipticity (1 - B/A) +real r #O Radius based on second moments +real cxx, cyy, cxy #O Ellipse parameters based on second moments +real aerr, berr, thetaerr #O Errors +real cxxerr, cyyerr, cxyerr #O Errors + +bool doerr +real x2, y2, xy, r2, d, f +real xvar, yvar, xycov, rvar, dvar, fvar + +begin + a = INDEFR + b = INDEFR + theta = INDEFR + elong = INDEFR + ellip = INDEFR + r = INDEFR + aerr = INDEFR + berr = INDEFR + thetaerr = INDEFR + cxxerr = INDEFR + cyyerr = INDEFR + cxyerr = INDEFR + + x2 = OBJ_X2(obj) + y2 = OBJ_Y2(obj) + xy = OBJ_XY(obj) + xvar = OBJ_XVAR(obj) + yvar = OBJ_YVAR(obj) + xycov = OBJ_XYCOV(obj) + + if (IS_INDEFR(x2) || IS_INDEFR(y2) || IS_INDEFR(xy)) + return + + r2 = x2 + y2 + if (r2 < 0.) + return + + doerr = !(IS_INDEF(xvar) || IS_INDEF(yvar) || IS_INDEF(xycov)) + if (doerr) { + rvar = xvar + yvar + if (rvar < 0.) + doerr = false + } + + r = sqrt (r2) + + d = x2 - y2 + theta = RADTODEG (atan2 (2 * xy, d) / 2.) + + if (doerr) { + dvar = xvar - yvar + thetaerr = atan2 (2 * xycov, dvar) / 2. + if (thetaerr < 0.) + thetaerr = INDEF + else + thetaerr = DEGTORAD (sqrt (thetaerr)) + } + + f = sqrt (d**2 + 4 * xy**2) + if (f > r2) + return + + if (doerr) { + fvar = sqrt (dvar**2 + 4 * xycov**2) + if (fvar > rvar) + doerr = false + } + + a = sqrt ((r2 + f) / 2) + b = sqrt ((r2 - f) / 2) + + if (doerr) { + aerr = sqrt ((rvar + fvar) / 2) + berr = sqrt ((rvar - fvar) / 2) + } + + ellip = 1 - b / a + if (b > 0.) + elong = a / b + + if (f == 0) { + cxx = 1. / (a * a) + cyy = 1. / (a * a) + cxy = 0 + } else { + cxx = y2 / f + cyy = x2 / f + cxy = -2 * xy / f + } + + if (doerr) { + if (fvar == 0) { + cxxerr = 1. / (aerr * aerr) + cyyerr = 1. / (berr * berr) + cxyerr = 0. + } else { + cxxerr = yvar / fvar + cyyerr = xvar / fvar + cxyerr = -2 * xycov / fvar + } + } + +end |