/* internal ana subroutines and functions */ /* mostly fft related stuff and solar mappings */ /*this is file fun5.c */ #include #include #include "ana_structures.h" #include #include "defs.h" /* for SGI only (?) */ #if __sgi #include #include #endif extern struct sym_desc sym[]; extern struct sym_list *subr_sym_list[]; extern int badmatch; extern int temp_base; extern int edb_context; extern int ana_float(), ana_double(); extern char *find_sym_name(); extern int ana_type_size[]; extern int vfix_top, num_ana_subr, next_user_subr_num; extern int lastmin_sym, lastmax_sym, max_temp; extern float float_arg(); extern double double_arg(); int fftdp = 0, revolve_mode = 0; float solar_rotate_a = 13.27; float solar_rotate_b = -1.68; float solar_rotate_c = -2.40; union types_ptr { byte *b; short *w; int *l; float *f; double *d;}; static int nold=0, *work, fftdpold, noldrdft=0, *workrdft, fftdprd; static int nolddcq, *workdcq, fftdpdcq; static double a1, a2, a3,a4,a5, a6, a7, a8, a9 ,a10, a11, a12, a13, a14, a15; static double subdx, subdy, xoff, yoff; void fluxcell(), subshift(); double sxvalue(), syvalue(), mert(), subshiftc(), subshiftc_apod(), meritc; /*------------------------------------------------------------------------- */ int ana_applyvfilter(narg,ps) /* routine for 3-D fft's */ /* specialized and complicated */ /* ana_applyvfilter,st,ct,vt,vs,vtop,vbot,blks,istart */ int narg, ps[]; { /* everybody must be defined already, ct and st must match the first dimensions of st and vt must match, the 2nd (and third if present) dimensions of st must match the dimensions of vs */ /* added a blks to allow only partial use of the ct and st array and istart, which together with blks, selects the relevant part of vs to use without doing an extraction before the call, these similify the call and improve performance at the price of specializing this routine. Also require st and ct to be 3-D. */ float *st, *ct, *vt, *vs, vtop, vbot, *pvt, *pvs, xq; int j, i, nd, nd2, nt, nt2, outer, outer2, n, iq, blks, istart, ny; struct ahead *h; /* setup (and check) pointers to the various floating point arrays */ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); if ( sym[iq].type != 3 ) return execute_error(120); h = (struct ahead *) sym[iq].spec.array.ptr; st = (float *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nt = h->dims[0]; if (nd != 3) { fprintf(stderr,"APPLYVFILTER - sine array must be 3-D\n"); return -1; } if (int_arg_stat(ps[6], &blks) != 1) return -1; if (int_arg_stat(ps[7], &istart) != 1) return -1; if (h->dims[2] < blks) { fprintf(stderr,"APPLYVFILTER - sine array 3rd dimension < blocking factor\n"); return -1; } outer = h->dims[1] * blks; iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(66); if ( sym[iq].type != 3 ) return execute_error(120); h = (struct ahead *) sym[iq].spec.array.ptr; ct = (float *) ((char *)h + sizeof(struct ahead)); nd2 = h->ndim; nt2 = h->dims[0]; if (nd2 != 3) { fprintf(stderr,"APPLYVFILTER - cosine array must be 3-D\n"); return -1; } if (nd2 > 1) for(j=1;jdims[j];} if (h->dims[2] < blks) { fprintf(stderr,"APPLYVFILTER - cosine array 3rd dimension < blocking factor\n"); return -1; } outer2 = h->dims[1] * blks; /* this routine is pretty specialized but still make consistency checks first ensure ct and st are structually the same */ /* the inner dimensions must agree and the product of the outers */ if (nt != nt2 || outer != outer2) { fprintf(stderr,"APPLYVFILTER - sine and cosine arrays don't match\n"); return -1; } iq = ps[2]; if ( sym[iq].class != 4 ) return execute_error(66); if ( sym[iq].type != 3 ) return execute_error(120); h = (struct ahead *) sym[iq].spec.array.ptr; vt = (float *) ((char *)h + sizeof(struct ahead)); nd2 = h->ndim; nt2 = h->dims[0]; if (nd2 > 1 || nt != nt2) { fprintf(stderr,"APPLYVFILTER - VT array doesn't match\n"); return -1; } iq = ps[3]; if ( sym[iq].class != 4 ) return execute_error(66); if ( sym[iq].type != 3 ) return execute_error(120); h = (struct ahead *) sym[iq].spec.array.ptr; vs = (float *) ((char *)h + sizeof(struct ahead)); /* here we use the entire first dimension and part of the second as specified by istart and blks, the result must = outers of st */ nd2 = h->ndim; if (nd2 != 2) { fprintf(stderr,"APPLYVFILTER - VS array must be 2-D\n"); return -1; } outer2 = h->dims[0]; if (h->dims[1] < (blks + istart)) { fprintf(stderr,"APPLYVFILTER - VS array, 2nd dimension < blocking range\n"); return -1; } outer2 *= blks; if ( outer != outer2) { printf("APPLYVFILTER - VS array doesn't match\n"); return -1; } /* now point to the part of vs array of interest */ vs += (h->dims[0]) * istart; /* now get the velocity limits */ if (float_arg_stat(ps[4], &vtop) != 1) return -1; if (float_arg_stat(ps[5], &vbot) != 1) return -1; /* note that we already got the last 2 arguments near the start */ while (outer--) { /* outer loop */ n = nt; pvs = vs++; pvt = vt; while (n--) { xq = (*pvt++) * (*pvs); if ( xq > vtop || xq < vbot ) { *ct = 0.0; *st = 0.0; } st++; ct++; } } return 1; } /*------------------------------------------------------------------------- */ int ana_fftwsc(narg,ps) /* sc routine */ /* trying the fftw libraries */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nn; int *wsave1, *wsave2; int ifac[20]; union types_ptr q1,q2,q3; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) iq = ana_float(1,ps); else iq = ana_double(1,ps); h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (nold != n || fftdp != fftdpold) { /* have to set up work space */ if (nold) free(work); /* delete old if any */ //mq = n * 28 +120; if ( fftdp == 0 ) mq = mq/2; //work = (int *) malloc( mq ); mq = 10 * n; if ( fftdp == 1 ) mq = mq * 2; wsave2 = (float *) malloc( mq ); #if NeXT if ( fftdp == 0 ) ezffti(&n, work); else ezfftid(&n, work); #else if ( fftdp == 0 ) ezffti_(&n, work); else ezfftid_(&n, work); // new_ezffti(n, wsave2, ifac); #endif nold = n; fftdpold = fftdp; } /* configure the output arrays */ nn = n/2 + 1; dim[0] = nn; iq = ps[1]; jq = ps[2]; if ( fftdp == 0 ) { type = 3; if (redef_array(iq, 3, nd, dim) != 1 || redef_array(jq, 3, nd, dim) != 1) return -1; } else { type = 4; if (redef_array(iq, 4, nd, dim) != 1 || redef_array(jq, 4, nd, dim) != 1) return -1; } /* get addresses */ h = (struct ahead *) sym[iq].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; q3.l = (int *) ((char *)h + sizeof(struct ahead)); /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT ezsc(&n, q1.f, q2.f, q3.f, work); #else ezsc_(&n, q1.f, q2.f, q3.f, work); #endif q1.f += nx; q2.f += nn; q3.f += nn; break; case 4: #if NeXT ezscd(&n, q1.d, q2.d, q3.d, work); #else ezscd_(&n, q1.d, q2.d, q3.d, work); #endif q1.d += nx; q2.d += nn; q3.d += nn; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_sc(narg,ps) /* sc routine */ /* sine-cosine style FFT -- sc, x, s, c */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nn; union types_ptr q1,q2,q3; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) iq = ana_float(1,ps); else iq = ana_double(1,ps); h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (nold != n || fftdp != fftdpold) { /* have to set up work space */ if (nold) free(work); /* delete old if any */ mq = n * 28 +120; if ( fftdp == 0 ) mq = mq/2; work = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) ezffti(&n, work); else ezfftid(&n, work); #else if ( fftdp == 0 ) ezffti_(&n, work); else ezfftid_(&n, work); #endif nold = n; fftdpold = fftdp; } /* configure the output arrays */ nn = n/2 + 1; dim[0] = nn; iq = ps[1]; jq = ps[2]; if ( fftdp == 0 ) { type = 3; if (redef_array(iq, 3, nd, dim) != 1 || redef_array(jq, 3, nd, dim) != 1) return -1; } else { type = 4; if (redef_array(iq, 4, nd, dim) != 1 || redef_array(jq, 4, nd, dim) != 1) return -1; } /* get addresses */ h = (struct ahead *) sym[iq].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; q3.l = (int *) ((char *)h + sizeof(struct ahead)); /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT ezsc(&n, q1.f, q2.f, q3.f, work); #else ezsc_(&n, q1.f, q2.f, q3.f, work); #endif q1.f += nx; q2.f += nn; q3.f += nn; break; case 4: #if NeXT ezscd(&n, q1.d, q2.d, q3.d, work); #else ezscd_(&n, q1.d, q2.d, q3.d, work); #endif q1.d += nx; q2.d += nn; q3.d += nn; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_power(narg,ps) /* power function */ /* power spectrum(s) */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, result_sym, nn; union types_ptr q1,q2,q3; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) iq = ana_float(1,ps); else iq = ana_double(1,ps); h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (nold != n || fftdp != fftdpold) { /* have to set up work space */ if (nold) free(work); /* delete old if any */ mq = n * 28 +120; if ( fftdp == 0 ) mq = mq/2; work = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) ezffti(&n, work); else ezfftid(&n, work); #else if ( fftdp == 0 ) ezffti_(&n, work); else ezfftid_(&n, work); #endif nold = n; fftdpold = fftdp; } /* configure the output array */ nn = n/2 + 1; dim[0] = nn; if ( fftdp == 0 ) { type = 3; result_sym = array_scratch(3, nd, dim); /*for the result */ } else { type = 4; result_sym = array_scratch(4, nd, dim); /*for the result */ } /* get address */ h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT ezpower(&n, q1.f, q2.f, work); #else ezpower_(&n, q1.f, q2.f, work); #endif q1.f += nx; q2.f += nn; break; case 4: #if NeXT ezpowerd(&n, q1.d, q2.d, work); #else ezpowerd_(&n, q1.d, q2.d, work); #endif q1.d += nx; q2.d += nn; break; } } return result_sym; } /*------------------------------------------------------------------------- */ int ana_scb(narg,ps) /* scb routine */ /* backwards sine-cosine style FFT -- scb, x, s, c */ int narg, ps[]; { static int nold, *work, fftdpold; int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nd2, outer2, nx2; union types_ptr q1,q2,q3; struct ahead *h; iq = ps[1]; jq = ps[2]; /* s and c must be arrays */ if ( sym[iq].class != 4 || sym[jq].class != 4) return execute_error(66); if ( fftdp == 0 ) { iq = ana_float(1, &iq); jq = ana_float(1, &jq); } else { iq = ana_double(1, &iq); jq = ana_double(1, &jq); } /* s and c must have similar structures */ h = (struct ahead *) sym[iq].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nx = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} h = (struct ahead *) sym[jq].spec.array.ptr; q3.l = (int *) ((char *)h + sizeof(struct ahead)); nd2 = h->ndim; nx2 = h->dims[0]; outer2 = 1; if (nd2 > 1) for(j=1;jdims[j];} /* the inner dimensions must agree and the product of the outers */ if (nx != nx2 || outer != outer2) { printf("SCB - sine and cosine arrays don't match\n"); return -1; } n = (nx - 1) * 2; /* scratch storage is needed, can we use the last one setup ? */ if (nold != n || fftdp != fftdpold) { /* have to set up work space */ if (nold) free(work); /* delete old if any */ mq = n * 28 +120; if ( fftdp == 0 ) mq = mq/2; work = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) ezffti(&n, work); else ezfftid(&n, work); #else if ( fftdp == 0 ) ezffti_(&n, work); else ezfftid_(&n, work); #endif nold = n; fftdpold = fftdp; } /* configure the output array */ dim[0] = n; iq = ps[0]; if ( fftdp == 0 ) { type = 3; if (redef_array(iq, 3, nd, dim) != 1) return -1; } else { type = 4; if (redef_array(iq, 4, nd, dim) != 1) return -1; } /* get addresse */ h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT ezscb(&n, q1.f, q3.f, q2.f, work); #else ezscb_(&n, q1.f, q3.f, q2.f, work); #endif q1.f += n; q2.f += nx; q3.f += nx; break; case 4: #if NeXT ezscbd(&n, q1.d, q3.d, q2.d, work); #else ezscbd_(&n, q1.d, q3.d, q2.d, work); #endif q1.d += n; q2.d += nx; q3.d += nx; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_rdft(narg,ps) /* discrete real transform routine */ /* real valued FFT -- rdft, x */ /* in place transform, maintain separate work areas from sc */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nn; union types_ptr q1,q2,q3; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) /* a replace is required since we need to output the same symbol */ { if (sym[iq].type != 3) ana_replace(iq, ana_float(1,&iq) ); } else { if (sym[iq].type != 4) ana_replace(iq, ana_double(1,&iq) ); } h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nx < 2) { printf("RDFT - vector length %d too short\n", nx); return -1;} if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (noldrdft != n || fftdp != fftdprd) { /* have to set up work space */ if (noldrdft) free(workrdft); /* delete old if any */ /* note, different size than sc and scb routines */ mq = n * 20 +120; if ( fftdp == 0 ) mq = mq/2; workrdft = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) rffti(&n, workrdft); else rfftid(&n, workrdft); #else if ( fftdp == 0 ) rffti_(&n, workrdft); else rfftid_(&n, workrdft); #endif noldrdft = n; fftdprd = fftdp; } if ( fftdp == 0 ) type = 3; else type = 4; /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT rfftf(&n, q1.f, workrdft); #else rfftf_(&n, q1.f, workrdft); #endif q1.f += nx; break; case 4: #if NeXT rfftfd(&n, q1.f, workrdft); #else rfftfd_(&n, q1.f, workrdft); #endif q1.d += nx; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_rdftb(narg,ps) /* discrete real transform reverse routine */ /* real valued back FFT -- rdftb, x */ /* in place transform, maintain separate work areas from sc */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nn; union types_ptr q1,q2,q3; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) /* a replace is required since we need to output the same symbol */ { if (sym[iq].type != 3) ana_replace(iq, ana_float(1,&iq) ); } else { if (sym[iq].type != 4) ana_replace(iq, ana_double(1,&iq) ); } h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nx < 4) { printf("RDFT - vector length %d too short\n", nx); return -1;} if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (noldrdft != n || fftdp != fftdprd) { /* have to set up work space */ if (noldrdft) free(workrdft); /* delete old if any */ /* note, different size than sc and scb routines */ mq = n * 20 +120; if ( fftdp == 0 ) mq = mq/2; workrdft = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) rffti(&n, workrdft); else rfftid(&n, workrdft); #else if ( fftdp == 0 ) rffti_(&n, workrdft); else rfftid_(&n, workrdft); #endif noldrdft = n; fftdprd = fftdp; } if ( fftdp == 0 ) type = 3; else type = 4; /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT rfftb(&n, q1.f, workrdft); #else rfftb_(&n, q1.f, workrdft); #endif q1.f += nx; break; case 4: #if NeXT rfftbd(&n, q1.f, workrdft); #else rfftbd_(&n, q1.f, workrdft); #endif q1.d += nx; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_dcqf(narg,ps) /* discrete cosine transform routine, f */ /* real valued cosine transform -- dcqf, x */ /* in place transform, maintain separate work areas from sc */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nn; union types_ptr q1; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) /* a replace is required since we need to output the same symbol */ { if (sym[iq].type != 3) ana_replace(iq, ana_float(1,&iq) ); } else { if (sym[iq].type != 4) ana_replace(iq, ana_double(1,&iq) ); } h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nx < 2) { printf("DCQF - vector length %d too short\n", nx); return -1;} if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (nolddcq != n || fftdp != fftdpdcq) { /* have to set up work space */ printf("setting up new workspace for n = %d\n", n); if (nolddcq) free(workdcq); /* delete old if any */ /* note, different size than sc and scb routines */ mq = n * 28 +120; if ( fftdp == 0 ) mq = mq/2; workdcq = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) cosqi(&n, workdcq); else cosqid(&n, workdcq); #else if ( fftdp == 0 ) cosqi_(&n, workdcq); else cosqid_(&n, workdcq); #endif nolddcq = n; fftdpdcq = fftdp; } if ( fftdp == 0 ) type = 3; else type = 4; /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT cosqf(&n, q1.f, workdcq); #else cosqf_(&n, q1.f, workdcq); #endif q1.f += nx; break; case 4: #if NeXT cosqfd(&n, q1.f, workdcq); #else cosqfd_(&n, q1.f, workdcq); #endif q1.d += nx; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_dcqb(narg,ps) /* discrete cosine transform routine, b */ /* real valued cosine transform -- dcqb, x */ /* in place transform, maintain separate work areas from sc */ int narg, ps[]; { int iq, nd, outer, nx, n, mq, dim[8], type, j, jq, nn; union types_ptr q1; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float the input */ if ( fftdp == 0 ) /* a replace is required since we need to output the same symbol */ { if (sym[iq].type != 3) ana_replace(iq, ana_float(1,&iq) ); } else { if (sym[iq].type != 4) ana_replace(iq, ana_double(1,&iq) ); } h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); /* find inner and product of all outer dimensions */ nd = h->ndim; nx = h->dims[0]; outer = 1; if (nx < 2) { printf("DCQB - vector length %d too short\n", nx); return -1;} if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j];} /* if nx is odd, we have to use nx-1 */ n = nx; if (n%2 != 0) n--; /* scratch storage is needed, can we use the last one setup ? */ if (nolddcq != n || fftdp != fftdpdcq) { /* have to set up work space */ printf("setting up new workspace for n = %d\n", n); if (nolddcq) free(workdcq); /* delete old if any */ /* note, different size than sc and scb routines */ mq = n * 28 +120; if ( fftdp == 0 ) mq = mq/2; workdcq = (int *) malloc( mq ); #if NeXT if ( fftdp == 0 ) cosqi(&n, workdcq); else cosqid(&n, workdcq); #else if ( fftdp == 0 ) cosqi_(&n, workdcq); else cosqid_(&n, workdcq); #endif nolddcq = n; fftdpdcq = fftdp; } if ( fftdp == 0 ) type = 3; else type = 4; /* loop over the outer dimensions */ while (outer--) { switch (type) { case 3: #if NeXT cosqb(&n, q1.f, workdcq); #else cosqb_(&n, q1.f, workdcq); #endif q1.f += nx; break; case 4: #if NeXT cosqbd(&n, q1.f, workdcq); #else cosqbd_(&n, q1.f, workdcq); #endif q1.d += nx; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_revolve(narg,ps) /* rotate a solar image */ int narg, ps[]; /* a subroutine - revolve,AIN,IXC,IYC,R,B,DT,xx,yy xx and yy are coordinates in original to use for revolved image */ { int iq, nd, nx, ny, result_sym, dim[8], jq; float r, b, dt, *ain, *xx, *yy, xc, yc; struct ahead *h; /* setup (and check) pointers to the various floating point arrays */ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; nd = h->ndim; nx = h->dims[0]; if (nd != 2) return execute_error(101); ny = h->dims[1]; dim[0] = nx; dim[1] = ny; //printf("nx, ny = %d %d\n", nx, ny); if (float_arg_stat(ps[1], &xc) != 1) return -1; if (float_arg_stat(ps[2], &yc) != 1) return -1; if (float_arg_stat(ps[3], &r) != 1) return -1; if (float_arg_stat(ps[4], &b) != 1) return -1; if (float_arg_stat(ps[5], &dt) != 1) return -1; iq = ps[6]; jq = ps[7]; nd = 2; dim[0] = nx; dim[1] = ny; if (redef_array(iq, 3, nd, dim) != 1 || redef_array(jq, 3, nd, dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; xx = (float *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; yy = (float *) ((char *)h + sizeof(struct ahead)); iq = revolve(nx, ny, xc, yc, nx, ny, xc, yc, r, b, dt, xx, yy); return 1; } /*------------------------------------------------------------------------- */ int revolve(nxin,nyin,xcin,ycin,nx,ny,ixc,iyc,r,b,dt,xx,yy) /* takes a solar image as seen from the earth and returns the location of each pixel (in arrays xx and yy) where it would be after the pixels rotated on the sun for a time dt; these can be used to construct an image of the sun after such a rotation; the size of the input image "ain" is [nxin,nyin], the pixel spacing is assumed uniform and the radius "r" of the sun is given in units of the pixel spacing in ain; e.g, if ain has an image point every 0.5", then the radius would be something like 2*960. The location of the center in the ain image is given by (xcin,ycin), also in units of the pixel spacing. These are cartesian coordinates from the (0,0) corner in the image with north and west positive. Frequently the solar center is outside the image of course. Inputting r and (xcin,ycin) this way establishes the scales with a minimum of input parameters. The size of the output image coordinates is given by [nx, ny]. The pixel scale is the same as ain, there is no current provision to change it. The desired disk position of the center of this image is given by (ixc,iyc), again in pixel spacing units and cartesian from disk center. r = radius in array index units b = solar b angle in degrees note that center coord. convention is the position of the center in array coord. base 1 (really not true, left over from fortran version?) centers outside the array are allowed the call from ana converts the user's center to this convention by adding 1 */ float *xx, *yy; int nyin, nxin, nx, ny; float xcin, ycin, ixc, iyc; float r, b, dt; { float rqd=1./57.2957795, sb, cb, dtr, xq, yq, rr, yq2, rq, cr, sc, slat, slon; float clat, rlon, clon, xp, yp, *pf1, *pf2; float latitude, th, dlon, sth; float sra, srb, src; int imin, jmin, imax, jmax, i, j, iq; sra = solar_rotate_a; srb = solar_rotate_b; src = solar_rotate_c; //printf("solar rotation values %f %f %f\n", sra, srb, src); //printf("inputs - nxin, nyin, xcin, ycin, nx, ny, ixc, iyc = %d %d %f %f %d %d %f %f\n", nxin, nyin, xcin, ycin, nx, ny, ixc, iyc); sb=sin(b*rqd); cb=cos(b*rqd); dtr=dt*rqd; /* always zero xx and yy first */ bzero( (char *) xx, nx*ny*sizeof(float)); bzero( (char *) yy, nx*ny*sizeof(float)); /* find column and row limits */ //printf("ixc, iyc, r = %f %f %f\n", ixc,iyc,r); imin = MAX( 0, (int) rint(ixc - r) ); jmin = MAX( 0, (int) rint(iyc - r) ); imax = MIN( nx, (int) rint(ixc+r) ); jmax = MIN( ny, (int) rint(iyc+r) ); rr=r*r; //printf("imin,imax,jmin,jmax = %d %d %d %d\n", imin,imax,jmin,jmax); /* the column loop */ for (j = jmin; j < jmax;j++) { yq = j - iyc; yq2 = yq * yq; /* note that yq is cos(chi)*r */ /* the row loop */ pf1 = xx + j*nx + imin; pf2 = yy + j*nx + imin; for (i = imin; i < imax;i++) { xq = i - ixc; rq = xq * xq + yq2; /* are we in the circle ? */ if (rq < rr) { /* step 1 - get lat. and long., compute sin(rho) and cos (rho) */ cr = sqrt(rr-rq); /* get sin(chi)*r */ sc = xq; slat = cr*sb+cb*yq; /* this is sin(lat.)*r */ xq = slat*slat; /* re-use xq variable */ clat = sqrt(rr-xq); /* and cos(lat.)*r */ if (ABS(clat) > 1.e-5) slon=sc/clat; else slon=0.0; /* slon is sin(long.) */ /* done with step 1, step 2 is to rotate back, include diff. rot. */ /* note that latitude doesn't change */ xq = xq/rr; /* now this is just slat*slat */ slon = MIN( slon, 1.0); slon = MAX( slon , -1.0); rlon = asin(slon); /* were we over a pole ? */ if (cr < (slat*sb)) rlon = 3.1415926 - rlon; if (revolve_mode == 0) { rlon = rlon - dtr * (sra + srb*xq + src*xq*xq); } else { /* only other option now is a special fit for Neal */ /* special fit for Neal, need actual latitude */ #if NeXT | defined(SOLARIS) | __APPLE__ latitude = asin(slat/r); #else latitude = asinf(slat/r); #endif th = 90. - ABS(latitude)*57.2957795; dlon =24.9848-1.19558*th+1.23576*th*th-0.0238393*th*th*th+0.000171238*th*th*th*th; /* and divide by sin(th) to get the angular rate */ #if NeXT | defined(SOLARIS) | __APPLE__ sth = sin(th*rqd); #else sth = sinf(th*rqd); #endif if (sth > 1.e-3) dlon=dlon/sth; else dlon = 0.0; /* that was m/s, convert to degrees per day and add to longitude */ rlon = rlon - dtr*.00711*dlon; } /* end of revolve_mode selection */ /* end of special fit for Neal */ /* end of step 2, step 3 is to get index in original */ /* check if not visible */ clon = cos(rlon); if ( (slat*sb+clat*cb*clon) > 0 ) { slon = sin(rlon); xp = clat * slon + xcin; /* go on only if this index is valid */ if ( (xp > 0) && (xp <= nxin)) { /* now do the column index */ yp = (slat * cb - clat * sb * clon) + ycin; if ( (yp > 0) && (yp <= nyin)) { *pf1 = xp; *pf2 = yp;} } } } pf1++; pf2++; } } return 1; } /*------------------------------------------------------------------------- */ int ana_mdi2earth(narg,ps) /* rotate a solar image */ int narg, ps[]; /* a subroutine - mdi2earth,AIN,xcmdi,ycmdi,rmdi,bmdi,clmdi,nx,ny,xc,yc,r,b,cl,xx,yy xx and yy are coordinates in original to use for remapped image note that AIN is used only to get the size of the original */ { int iq, nd, nx, ny, result_sym, dim[8], jq, nxmdi,nymdi; float xcmdi,ycmdi,rmdi,bmdi,clmdi,xc,yc,r,b,cl, *xx, *yy; struct ahead *h; /* setup (and check) pointers to the various floating point arrays */ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; nd = h->ndim; nxmdi = h->dims[0]; if (nd != 2) return execute_error(101); nymdi = h->dims[1]; printf("nxmdi, nymdi = %d %d\n", nxmdi, nymdi); if (float_arg_stat(ps[1], &xcmdi) != 1) return -1; if (float_arg_stat(ps[2], &ycmdi) != 1) return -1; if (float_arg_stat(ps[3], &rmdi) != 1) return -1; if (float_arg_stat(ps[4], &bmdi) != 1) return -1; if (float_arg_stat(ps[5], &clmdi) != 1) return -1; if (int_arg_stat(ps[6], &nx) != 1) return -1; if (int_arg_stat(ps[7], &ny) != 1) return -1; if (float_arg_stat(ps[8], &xc) != 1) return -1; if (float_arg_stat(ps[9], &yc) != 1) return -1; if (float_arg_stat(ps[10], &r) != 1) return -1; if (float_arg_stat(ps[11], &b) != 1) return -1; if (float_arg_stat(ps[12], &cl) != 1) return -1; iq = ps[13]; jq = ps[14]; dim[0] = nx; dim[1] = ny; nd = 2; dim[0] = nx; dim[1] = ny; printf("nx, ny = %d %d\n", nx, ny); if (redef_array(iq, 3, nd, dim) != 1 || redef_array(jq, 3, nd, dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; xx = (float *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; yy = (float *) ((char *)h + sizeof(struct ahead)); iq = mdi2earth(nxmdi,nymdi,xcmdi,ycmdi,rmdi,bmdi,clmdi,nx,ny,xc,yc,r,b,cl,xx,yy); return 1; } /*------------------------------------------------------------------------- */ int mdi2earth(nxmdi,nymdi,xcmdi,ycmdi,rmdi,bmdi,clmdi,nx,ny,xc,yc,r,b,cl,xx,yy) /* use to change the point of view for a solar observation */ float *xx, *yy; int nymdi, nxmdi, nx, ny; float xcmdi, ycmdi, xc, yc; float r, rmdi, b, bmdi, cl, clmdi; { float rqd=1./57.2957795, sb, cb, xq, yq, rr, yq2, rq, cr, sc, slat, slon; float clat, rlon, clon, xp, yp, *pf1, *pf2; float sbmdi, cbmdi, delcl, r_rat; int imin, jmin, imax, jmax, i, j, iq; /* xc, yc, r, b, cl are for the output image */ sb = sin(b*rqd); cb = cos(b*rqd); sbmdi = sin(bmdi*rqd); cbmdi = cos(bmdi*rqd); delcl = (cl - clmdi)*rqd; r_rat = rmdi/r; /* xx and yy are dimensioned for the output image */ /* always zero xx and yy first */ bzero( (char *) xx, nx*ny*sizeof(float)); bzero( (char *) yy, nx*ny*sizeof(float)); /* find column and row limits based on a box enclosing the disk, we check for the "corners" inside the loop, this reduces the work if the output space is larger than the sun */ imin = MAX( 0, (int) rint(xc - r) ); jmin = MAX( 0, (int) rint(yc - r) ); imax = MIN( nx, (int) rint(xc+r) ); jmax = MIN( ny, (int) rint(yc+r) ); rr=r*r; /* the column loop */ for (j = jmin; j < jmax;j++) { yq = j - yc; yq2 = yq * yq; /* note that yq is cos(chi)*r */ /* the row loop */ pf1 = xx + j*nx + imin; pf2 = yy + j*nx + imin; for (i = imin; i < imax;i++) { xq = i - xc; rq = xq * xq + yq2; /* are we in the circle ? */ if (rq < rr) { /* step 1 - get lat. and long for the output image pixels */ cr = sqrt(rr-rq); slat = cr*sb+cb*yq; /* this is sin(lat.)*r */ clat = sqrt(rr-slat*slat); /* and cos(lat.)*r */ /* the longitude is asin( x/(r*cos(lat))) */ if (ABS(clat) > 1.e-5) slon=xq/clat; else slon=0.0; /* slon is sin(long.) */ slon = MIN( slon, 1.0); slon = MAX( slon , -1.0); rlon = asin(slon); /* were we over a pole ? */ if (cr < (slat*sb)) rlon = 3.1415926 - rlon; /* we now effectively have the latitude and Carrington longitude for the output image for this pixel position. Actually we have the r*sin(lat), r*cos(lat) and CL - CLl0 (in rlon). */ rlon = rlon + delcl; /* rlon is now the longitude wrt y axis for the MDI image */ clon = cos(rlon); slon = sin(rlon); /* to get coordinates in the MDI image, we have to apply the radius seen by MDI and use the MDI B angle, the CL correction already done just above */ xp = r_rat * clat * slon + xcmdi; /* note that x doesn't depend on B here */ /* go on only if this index is valid */ if ( (xp > 0) && (xp <= nxmdi)) { /* now do the column index */ yp = r_rat * (slat * cbmdi - clat * sbmdi * clon) + ycmdi; if ( (yp > 0) && (yp <= nymdi)) { *pf1 = xp; *pf2 = yp;} } } pf1++; pf2++; } } return 1; } /*------------------------------------------------------------------------- */ int ana_postel(narg,ps) /* return projection coordinates */ int narg, ps[]; /* postel, nx, ny, scale, r, b, phi, latc, lonc, xc, yc, xx, yy see sunpostel below for definitions */ { int iq, jq, nd, nx, ny, dim[8]; float r, b, scale, phi, latc, lonc, xc, yc, *xx, *yy; struct ahead *h; if (int_arg_stat(ps[0], &nx) != 1) return -1; if (int_arg_stat(ps[1], &ny) != 1) return -1; if (float_arg_stat(ps[2], &scale) != 1) return -1; if (float_arg_stat(ps[3], &r) != 1) return -1; if (float_arg_stat(ps[4], &b) != 1) return -1; if (float_arg_stat(ps[5], &phi) != 1) return -1; if (float_arg_stat(ps[6], &latc) != 1) return -1; if (float_arg_stat(ps[7], &lonc) != 1) return -1; if (float_arg_stat(ps[8], &xc) != 1) return -1; if (float_arg_stat(ps[9], &yc) != 1) return -1; iq = ps[10]; jq = ps[11]; nd = 2; dim[0] = nx; dim[1] = ny; if (redef_array(iq, 3, nd, dim) != 1 || redef_array(jq, 3, nd, dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; xx = (float *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; yy = (float *) ((char *)h + sizeof(struct ahead)); sunpostel(nx, ny, scale, r, b, phi, latc, lonc, xc, yc, xx, yy); } /*------------------------------------------------------------------------- */ int sunpostel(nx, ny, scale, r, b, phi, latc, lonc, xc, yc, xx, yy) /* maps a solar image via a Postel projection, the pixel grid for the final projection is given by (nx, ny) and scale (assumed the same in both directions) in units of degrees; for small ranges of angles, the result has nearly square pixels in angle units. The B angle b and a roll angle phi are also input; also the central longitude of the area to be projected, this effectively rotates the image by this amount and then projects the area at disk center. A central latitude is also input to allow an area not centered at the equator. Together, this means that the point at (latc,longc) on the original image gets mapped to the center of the projection. Also input is the radius of the sun in units of the original image pixels and the location of the center in the same units. Returned are arrays of the x and y indices for the original image that map to the projection. */ int nx, ny; float scale, r, b, phi, lonc, latc, xc, yc, *xx, *yy; { float rqd=1./57.2957795, sb, cb, dtr, xq, yq, rr, yq2, rq, cr, sc, slat, slon; float clat, rlon, clon, lat, lon, p, cphi, sphi, sinr,sinp, lq, aq, bq; float slatc, clatc; int i, j, iq, nx2, ny2; xq = b * rqd; sb=sin(xq); cb=cos(xq); slatc = sin(latc * rqd); clatc = cos(latc * rqd); scale = scale * rqd; /* convert scale to radians */ xq = phi * rqd; cphi = cos(xq) * r; sphi = sin(xq) * r; nx2 = nx/2; ny2 = ny/2; /* the column loop */ for (j = 0; j < ny; j++) { yq = (float) (j - ny2); yq2 = yq * yq; /* the row loop */ for (i = 0; i < nx; i++) { xq = (float) (i - nx2); rq = scale* sqrt(xq * xq + yq2); p = atan2(yq, xq); if (xq == 0 && yq == 0) p = 0.0; /* step 1 - get lat. and long., postel projection */ sinr = sin(rq); sinp = sin(p); slat = sinr * sinp * clatc + cos(rq) * slatc; clat = sqrt(1. - slat*slat); lon = asin( sinr * cos(p)/ clat); /* were we over a pole ? */ if (cos(rq) < (slat*sb)) lon = 3.1415926 - lon; /* first get the coords in units of the radius without phi */ lq = lon + lonc * rqd; clon = cos(lq); slon = sin(lq); /* note that x is independent of B if phi is 0 */ aq = clat * slon; bq = slat * cb - clat *clon * sb; /* and rotate by phi, include r, and the center offset */ *xx++ = aq * cphi - bq * sphi + xc; *yy++ = bq * cphi + aq * sphi + yc; } } return 1; } /*------------------------------------------------------------------------- */ int ana_latlongrid(narg,ps) /* return coordinates for lat/lon grid*/ int narg, ps[]; /* postel, nx, ny, scale, r, b, phi, latc, lonc, xc, yc, xx, yy see sunlatlong below for definitions */ /* this is identical to ana_sunpostel except for the latc line, could be combined to shorten code */ { int iq, jq, nd, nx, ny, dim[8]; float r, b, scale, phi, latc, lonc, xc, yc, *xx, *yy; struct ahead *h; if (int_arg_stat(ps[0], &nx) != 1) return -1; if (int_arg_stat(ps[1], &ny) != 1) return -1; if (float_arg_stat(ps[2], &scale) != 1) return -1; if (float_arg_stat(ps[3], &r) != 1) return -1; if (float_arg_stat(ps[4], &b) != 1) return -1; if (float_arg_stat(ps[5], &phi) != 1) return -1; if (float_arg_stat(ps[6], &latc) != 1) return -1; if (float_arg_stat(ps[7], &lonc) != 1) return -1; if (float_arg_stat(ps[8], &xc) != 1) return -1; if (float_arg_stat(ps[9], &yc) != 1) return -1; iq = ps[10]; jq = ps[11]; nd = 2; dim[0] = nx; dim[1] = ny; if (redef_array(iq, 3, nd, dim) != 1 || redef_array(jq, 3, nd, dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; xx = (float *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; yy = (float *) ((char *)h + sizeof(struct ahead)); sunlatlong(nx, ny, scale, r, b, phi, latc, lonc, xc, yc, xx, yy); } /*------------------------------------------------------------------------- */ int sunlatlong(nx, ny, scale, r, b, phi, latc, lonc, xc, yc, xx, yy) /* maps a solar image onto a latitude/longitude grid, the pixel grid for the final projection is given by (nx, ny) and scale (assumed the same in both directions) in units of degrees; the result has square pixels in angle units. The B angle b and a roll angle phi are also input; also the central longitude of the area to be projected, this effectively rotates the image by this amount and then projects the area at disk center. A central latitude is also input to allow an area not centered at the equator. Together, this means that the point at (latc,longc) on the original image gets mapped to the center of the grid. Also input is the radius of the sun in units of the original image pixels and the location of the center in the same units. Returned are arrays of the x and y indices for the original image that correspond to the latitude/longitude grid. This was modified from the more complicated sunpostel routine */ int nx, ny; float scale, r, b, phi, lonc, latc, xc, yc, *xx, *yy; { float rqd=1./57.2957795, sb, cb, dtr, xq, yq, rr, yq2, rq, cr, sc, slat, slon; float clat, rlon, clon, ylat, xlon, p, cphi, sphi, sinr,sinp, lq, aq, bq; float slatc, clatc, nxc, nyc; int i, j, iq, nx2, ny2; xq = b * rqd; sb=sin(xq); cb=cos(xq); latc = latc * rqd; lonc = lonc * rqd; scale = scale * rqd; /* convert scale to radians */ xq = phi * rqd; cphi = cos(xq) * r; sphi = sin(xq) * r; nx2 = nx/2; ny2 = ny/2; nxc = (nx-1.0)*.5; nyc = (ny-1.0)*.5; /* the column loop */ for (j = 0; j < ny; j++) { ylat = (float) (j - nyc)*scale + latc; clat = cos(ylat); slat = sin(ylat); /* the row loop */ for (i = 0; i < nx; i++) { xlon = (float) (i - nxc)*scale + lonc; clon = cos(xlon); slon = sin(xlon); /* note that x is independent of B if phi is 0 */ aq = clat * slon; bq = slat * cb - clat *clon * sb; /* and rotate by phi, include r, and the center offset */ *xx++ = aq * cphi - bq * sphi + xc; *yy++ = bq * cphi + aq * sphi + yc; } } return 1; } /*------------------------------------------------------------------------- */ int ana_subshift(narg,ps) /* LCT for a cell */ /* wants 2 arrays, already F*8 and extracted, both the same size */ /* returns the shift */ int narg, ps[]; { int iq, jq; struct ahead *h; double *x1, *x2; int nx, ny; iq =ps[0]; jq =ps[1]; if (sym[iq].type != 4 || sym[jq].type != 4) { printf("ANA_SUBSHIFT: input arrays must be F*8\n"); return 0; } if (sym[iq].class != 4 || sym[jq].class != 4) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; x1 = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); nx = h->dims[0]; ny = h->dims[1]; h = (struct ahead *) sym[jq].spec.array.ptr; x2 = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); if ( nx != h->dims[0] || ny != h->dims[1]) return execute_error(103); subshift(x1, x2, nx, ny); /* expect result in xoff and yoff */ if (redef_scalar( ps[2], 4, &xoff) != 1) return execute_error(13); if (redef_scalar( ps[3], 4, &yoff) != 1) return execute_error(13); return 1; } /*------------------------------------------------------------------------- */ int ana_subshiftc(narg,ps) /* LCT for a cell, sym version */ /* subshiftc, s1b, s2b, xoff, yoff, mask */ /* 3/4/97 added apodizer array, optional */ /* wants 2 or 3 arrays, already F*8 and extracted, both the same size */ /* returns the shift */ int narg, ps[]; { int iq, jq, kq; struct ahead *h; double *x1, *x2, *msk; int nx, ny; iq =ps[0]; jq =ps[1]; if (sym[iq].type != 4 || sym[jq].type != 4 ) { printf("ANA_SUBSHIFT: input arrays must be F*8\n"); return 0; } if (sym[iq].class != 4 || sym[jq].class != 4 ) return execute_error(66); /* if a mask, also some rules */ if (narg > 4) { kq =ps[4]; if (sym[kq].type != 4) { printf("ANA_SUBSHIFT: input arrays must be F*8\n"); return 0; } if (sym[kq].class != 4) return execute_error(66); } h = (struct ahead *) sym[iq].spec.array.ptr; x1 = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); nx = h->dims[0]; ny = h->dims[1]; h = (struct ahead *) sym[jq].spec.array.ptr; x2 = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); if ( nx != h->dims[0] || ny != h->dims[1]) return execute_error(103); if (narg > 4) { h = (struct ahead *) sym[kq].spec.array.ptr; msk = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); if ( (nx-1) != h->dims[0] || (ny-1) != h->dims[1]) return execute_error(103); } switch (narg) { case 4: meritc = subshiftc(x1, x2, nx, ny); break; case 5: meritc = subshiftc_apod(x1, x2, msk, nx, ny); break; default: printf("internal error in ana_subshiftc\n"); return 0; } /* expect result in xoff and yoff */ if (redef_scalar( ps[2], 4, &xoff) != 1) return execute_error(13); if (redef_scalar( ps[3], 4, &yoff) != 1) return execute_error(13); return 1; } /*------------------------------------------------------------------------- */ int ana_fluxcell(narg,ps) /* optical flux for a cell */ /* wants 3 arrays, first 2 are the subimages and third is mask */ /* returns the shift */ int narg, ps[]; { int iq, jq, kq; struct ahead *h; double *x1, *x2, *msk; int nx, ny; iq =ps[0]; jq =ps[1]; kq =ps[2]; if (sym[iq].class != 4 || sym[jq].class != 4 || sym[kq].class != 4) return execute_error(66); /* always convert to dp, at least for now */ iq = ana_double(1,&ps[0]); jq = ana_double(1,&ps[1]); kq = ana_double(1,&ps[2]); h = (struct ahead *) sym[iq].spec.array.ptr; x1 = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); nx = h->dims[0]; ny = h->dims[1]; h = (struct ahead *) sym[jq].spec.array.ptr; x2 = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); if ( nx != h->dims[0] || ny != h->dims[1]) return execute_error(103); h = (struct ahead *) sym[kq].spec.array.ptr; msk = (double *) ((char *)h + sizeof(struct ahead)); if ( h->ndim != 2 ) return execute_error(101); if ( nx != h->dims[0] || ny != h->dims[1]) return execute_error(103); /* got pointers to everybody, now let fluxcell do the work */ fluxcell(x1, x2, msk, nx, ny); /* expect result in xoff and yoff */ if (redef_scalar( ps[3], 4, &xoff) != 1) return execute_error(13); if (redef_scalar( ps[4], 4, &yoff) != 1) return execute_error(13); return 1; } /*------------------------------------------------------------------------- */ void fluxcell(xa, xb, msk, nx, ny) double *xa, *xb, *msk; int nx, ny; { /* unfinished */ double *x1, *x2, *x3, *x4; double dx, dt, dx2; int j; x1 = xa + nx; x2 = xb + nx; /* outer j loop of ny-2 length */ j = ny - 2; while (j--) { dx = *(x1+2) - *x1; dx2 = *(x2+2) - *x2; dt = *(x2+1) - *(x1+1); x1++; x2++; } } /*------------------------------------------------------------------------- */ double mert(sx,sy) double sx,sy; { double w0, w1,w2,w3,mq; w0 = (1-sx)*(1-sy); w1 = sx*(1-sy); w2 = sy*(1-sx); w3 = sx*sy; mq = w0*w0*a1 + w1*w1*a2 + w2*w2*a3 + w3*w3*a4 + w0*w1*a5 + w0*w2*a6 + w0*w3*a7 + w1*w2*a8 + w1*w3*a9 + w2*w3*a10 + w0*a11 + w1*a12 + w2*a13 + w3*a14 + a15; return mq; } /*------------------------------------------------------------------------- */ double mertc(sx,sy) double sx,sy; { double xq; xq = a1 + sx*sx*a2 + sy*sy*a3 + sx*sx*sy*sy*a4 + sx*a5 + sy*a6 + sx*sy*a7 + sx*sx*sy*a9 + sx*sy*sy*a10; return xq; } /*------------------------------------------------------------------------- */ double sxvalue(sy) double sy; { double syc, c1, c2, xq; syc = 1 - sy; c1 = syc*syc*(a5-2.0*a1) +syc*sy*(a7-2.0*a6+a8) +sy*sy*(a10-2.0*a3) +syc*(a12-a11) +sy*(a14-a13); c2 = syc*syc*(a1+a2-a5) +sy*sy*(a3+a4-a10) +syc*sy*(a6-a7-a8+a9); if (c2 == 0.0) { /* no longer print this, it happens for 0 counts, instead bump up badmatch */ /* printf("sxvalue singular, c1, c2 = %g %g\n", c1,c2);*/ badmatch++; return -1.0;} if ( isnan(c1) == 1 || isnan(c2) ==1 ) { printf("sxvalue problem:\n"); printf("sy, syc = %14.6e %14.6e\n", sy,syc); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10,a11 = %14.6e %14.6e %14.6e %14.6e\n",a8,a9,a10,a11); printf("a12,a13,a14,a15 = %14.6e %14.6e %14.6e %14.6e\n",a12,a13,a14,a15); return -1.0; } xq = -.5*c1/c2; return xq; } /*------------------------------------------------------------------------- */ double syvalue(sx) double sx; { double sxc, c1, c2, xq; sxc = 1 - sx; c1 = sxc*sxc*(a6-2.0*a1) +sx*sx*(a9-2.0*a2) +sxc*sx*(a7-2.0*a5+a8) +sxc*(a13-a11) +sx*(a14-a12); c2 = sxc*sxc*(a1+a3-a6) +sx*sx*(a2+a4-a9) +sxc*sx*(a5-a7-a8+a10); if (c2 == 0.0) { /* no longer print this, it happens for 0 counts, instead bump up badmatch */ /* printf("syvalue singular, c1, c2 = %g %g\n",c1,c2); */ badmatch++; return -1.0;} if ( isnan(c1) == 1 || isnan(c2) ==1 ) { printf("syvalue problem:\n"); printf("sx, sxc = %14.6e %14.6e\n", sx,sxc); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10,a11 = %14.6e %14.6e %14.6e %14.6e\n",a8,a9,a10,a11); printf("a12,a13,a14,a15 = %14.6e %14.6e %14.6e %14.6e\n",a12,a13,a14,a15); return -1.0; } xq = -.5*c1/c2; return xq; } /*------------------------------------------------------------------------- */ void getsxsy() { /* iterate to the min (if any), loading results in globals subdx and subdy */ /* seed with syz = 0.5 */ int n; double sxz, syz, delx, dely; syz=sxz=.5; n = 11; while (n--) { /* get new sxz */ sxz = sxvalue(syz); if ( sxz < -0.5 || sxz > 1.5) { subdx = subdy = 0.0; return; } /* and get a new syz */ syz = syvalue(sxz); if ( syz < -0.5 || syz > 1.5) { subdx = subdy = 0.0; return; } /* printf(" sxz, syz = %g %g\n", sxz,syz); */ } /* for small shifts, let values very close to the line but on the other side remain as contenders */ /* 4/30/96, try another scheme, let them have a value along the line and then we will compare on merit, this will allow 0's to be a result when the merit is really at a min along the line */ if (sxz >= 0.0 && sxz <=1.0 && syz >= 0.0 && syz <=1.0 ) { /* these OK */ subdx = sxz; subdy = syz; } else { /* maybe OK to use a zero or even a 1 if rather close to the line */ if (sxz >= -0.1 && sxz <=1.1 && syz >= -0.1 && syz <=1.1 ) { /* just range limit them, this should sort it out */ subdx = MAX(subdx, 0.0); subdy = MAX(subdy, 0.0); subdx = MIN(subdx, 1.0); subdy = MIN(subdy, 1.0); } else { subdx =0; subdy =0; } } /*printf("final sxz, syz = %g %g\n", sxz,syz);*/ return; } /*------------------------------------------------------------------------- */ void subshift(x, r, nx, ny) /* q is the reference cell (old m1), r is the one we interpolate (old m2) */ /* we assume that the arrays have already been converted to F*8 */ double *r, *x; int nx, ny; { int nxs, nys; double sum, sumx, rs1, rs0, rs2, rs3, cs0, cs1, cs2, cs3, t2, t1, t0, t3; double parts[5][5], xx[3][3], xdx[3][2], xdy[2][3], xmmpp[2][2], xppmm[2][2]; double partsdx[5][3], partsdy[3][5], partsppmm[3][3], partsmmpp[3][3]; double cmm,c0m,cpm,cm0,c00,cp0,cmp,c0p,cpp,sumxx; double qbest, qcur, outside, qd; int i, j, nxm2, nym2, ii, jj, mflag; double *rp, *rp2, *row, *rp3, *rp0, *rowq, *rp1, *qp; nxs = nx; nys = ny; nxm2 = nx - 2; nym2 = ny - 2; /* the sums of squares */ /* first rearrange as 25 pieces in an array of 5x5, mostly to have it semi organized, may be faster to do it differently */ rp = r; t0 = *rp++; parts[0][0] = t0*t0; t0 = *rp++; parts[0][1] = t0*t0; i = nx-4; sum = 0.0; while (i--) { t0 = *rp++; sum += t0 * t0; } parts[0][2] = sum; t0 = *rp++; parts[0][3] = t0*t0; t0 = *rp++; parts[0][4] = t0*t0; rp = r+nxs; t0 = *rp++; parts[1][0] = t0*t0; t0 = *rp++; parts[1][1] = t0*t0; i = nx-4; sum = 0.0; while (i--) { t0 = *rp++; sum += t0 * t0; } parts[1][2] = sum; t0 = *rp++; parts[1][3] = t0*t0; t0 = *rp++; parts[1][4] = t0*t0; /* third row has column sums and center sum */ row = r + nxs+ nxs; j = ny-4; cs0 = cs1 = cs2 = cs3 = sum = 0.0; while (j--) { rp = row; t0 = *rp++; row = row + nxs; cs0 += t0 * t0; t0 = *rp++; cs1 += t0 * t0; i = nx-4; while (i--) { t0 = *rp++; sum += t0 * t0; } t0 = *rp++; cs2 += t0 * t0; t0 = *rp++; cs3 += t0 * t0; } parts[2][0] = cs0; parts[2][1] = cs1; parts[2][2] = sum; parts[2][3] = cs2; parts[2][4] = cs3; rp = row; t0 = *rp++; parts[3][0] = t0*t0; t0 = *rp++; parts[3][1] = t0*t0; i = nx-4; sum = 0.0; while (i--) { t0 = *rp++; sum += t0 * t0; } parts[3][2] = sum; t0 = *rp++; parts[3][3] = t0*t0; t0 = *rp++; parts[3][4] = t0*t0; rp = row+nxs; t0 = *rp++; parts[4][0] = t0*t0; t0 = *rp++; parts[4][1] = t0*t0; i = nx-4; sum = 0.0; while (i--) { t0 = *rp++; sum += t0 * t0; } parts[4][2] = sum; t0 = *rp++; parts[4][3] = t0*t0; t0 = *rp++; parts[4][4] = t0*t0; /* now each sum of squares is a sum of 9 parts */ for (j=0;j<3;j++) { for (i=0;i<3;i++) { xx[j][i] = 0.0; for (jj=0;jj<3;jj++) { for (ii=0;ii<3;ii++) { xx[j][i] += parts[j+jj][i+ii]; }}}} /* the cross term cases with delta 1 in the x direction */ /* we build 15 sums which are combined into 6 results */ rp = r; t0 = *rp++; t1 = *rp++; partsdx[0][0] = t0 * t1; i = nx-3; sum = 0.0; while (i--) { t0 = t1; t1 = *rp++; sum += t0 * t1; } partsdx[0][1] = sum; /* one more for the top row */ t0 = t1; t1 = *rp++; partsdx[0][2] = t0 * t1; /* second row much the same */ rp = r+nxs; t0 = *rp++; t1 = *rp++; partsdx[1][0] = t0 * t1; i = nx-3; sum = 0.0; while (i--) { t0 = t1; t1 = *rp++; sum += t0 * t1; } partsdx[1][1] = sum; t0 = t1; t1 = *rp++; partsdx[1][2] = t0 * t1; /* third row has larger areas, 2 columns and the center */ row = r + nxs + nxs; j = ny-4; cs0 = cs1 = sum = 0.0; while (j--) { rp = row; t0 = *rp++; t1 = *rp++; row = row + nxs; cs0 += t0 * t1; /* left column sum */ i = nx-3; while (i--) { t0 = t1; t1 = *rp++; sum += t0 * t1; } t0 = t1; t1 = *rp++; cs1 += t0 * t1; } partsdx[2][0] = cs0; partsdx[2][1] = sum; partsdx[2][2] = cs1; /* the last 2 rows are just one high again */ rp = row; t0 = *rp++; t1 = *rp++; partsdx[3][0] = t0 * t1; i = nx-3; sum = 0.0; while (i--) { t0 = t1; t1 = *rp++; sum += t0 * t1; } partsdx[3][1] = sum; /* one more for the top row */ t0 = t1; t1 = *rp++; partsdx[3][2] = t0 * t1; /* second row much the same */ rp = row+nxs; t0 = *rp++; t1 = *rp++; partsdx[4][0] = t0 * t1; i = nx-3; sum = 0.0; while (i--) { t0 = t1; t1 = *rp++; sum += t0 * t1; } partsdx[4][1] = sum; t0 = t1; t1 = *rp++; partsdx[4][2] = t0 * t1; /* now each sum of cross terms is a sum of 9 parts */ for (j=0;j<3;j++) { for (i=0;i<2;i++) { xdx[j][i] = 0.0; for (jj=0;jj<3;jj++) { for (ii=0;ii<2;ii++) { xdx[j][i] += partsdx[j+jj][i+ii]; }}}} /* the terms with a simple shift in y are similar in concept */ /* we build 15 sums which are combined into 6 results */ rp = r; rp2 = rp+nxs; partsdy[0][0] = *rp++ * *rp2++; partsdy[0][1] = *rp++ * *rp2++; i = nx-4; sum = 0.0; while (i--) { sum += *rp++ * *rp2++; } partsdy[0][2] = sum; partsdy[0][3] = *rp++ * *rp2++; partsdy[0][4] = *rp++ * *rp2++; /* second row has column sums and center sum */ row = r + nxs; j = ny-3; cs0 = cs1 = cs2 = cs3 = sum = 0.0; while (j--) { rp = row; rp2 = rp+nxs; row = row + nxs; cs0 += *rp++ * *rp2++; /* left column sum */ cs1 += *rp++ * *rp2++; /* second column sum */ i = nx-4; while (i--) { sum += *rp++ * *rp2++; } /* center */ cs2 += *rp++ * *rp2++; /* second last column sum */ cs3 += *rp++ * *rp2++; /* last column sum */ } partsdy[1][0] = cs0; partsdy[1][1] = cs1; partsdy[1][2] = sum; partsdy[1][3] = cs2; partsdy[1][4] = cs3; /* the last row */ rp = row; rp2 = rp+nxs; partsdy[2][0] = *rp++ * *rp2++; partsdy[2][1] = *rp++ * *rp2++; i = nx-4; sum = 0.0; while (i--) { sum += *rp++ * *rp2++; } partsdy[2][2] = sum; partsdy[2][3] = *rp++ * *rp2++; partsdy[2][4] = *rp++ * *rp2++; /* now each sum of cross terms is a sum of 9 parts */ for (j=0;j<2;j++) { for (i=0;i<3;i++) { xdy[j][i] = 0.0; for (jj=0;jj<2;jj++) { for (ii=0;ii<3;ii++) { xdy[j][i] += partsdy[j+jj][i+ii]; }}}} /* diagonal terms, 2 types, first the larger gap, mmpp type */ /* we build 9 sums which are combined into 4 results */ rp = r; rp2 = rp+nxs+1; partsmmpp[0][0] = *rp++ * *rp2++; i = nx-3; sum = 0.0; while (i--) { sum += *rp++ * *rp2++; } partsmmpp[0][1] = sum; partsmmpp[0][2] = *rp++ * *rp2++; row = r + nxs; j=ny-3; cs0 = cs1 = sum = 0.0; while (j--) { rp = row; rp2 = rp+nxs+1; row = row + nxs; cs0 += *rp++ * *rp2++; /* left column sum */ i = nx-3; while (i--) { sum += *rp++ * *rp2++; } cs1 += *rp++ * *rp2++; /* right column sum */ } partsmmpp[1][0] = cs0; partsmmpp[1][1] = sum; partsmmpp[1][2] = cs1; rp = row; rp2 = rp+nxs+1; partsmmpp[2][0] = *rp++ * *rp2++; i = nx-3; sum = 0.0; while (i--) { sum += *rp++ * *rp2++; } partsmmpp[2][1] = sum; partsmmpp[2][2] = *rp++ * *rp2++; /* now each sum of cross terms is a sum of 4 parts */ for (j=0;j<2;j++) { for (i=0;i<2;i++) { xmmpp[j][i] = 0.0; for (jj=0;jj<2;jj++) { for (ii=0;ii<2;ii++) { xmmpp[j][i] += partsmmpp[j+jj][i+ii]; }}}} /* diagonal terms, 2 types, first the larger gap, mmpp type */ /* we build 9 sums which are combined into 4 results */ rp = r+1; rp2 = rp+nxs-1; partsppmm[0][0] = *rp++ * *rp2++; i = nx-3; sum = 0.0; while (i--) { sum += *rp++ * *rp2++; } partsppmm[0][1] = sum; partsppmm[0][2] = *rp++ * *rp2++; row = r + nxs+1; j=ny-3; cs0 = cs1 = sum = 0.0; while (j--) { rp = row; rp2 = rp+nxs-1; row = row + nxs; cs0 += *rp++ * *rp2++; /* left column sum */ i = nx-3; while (i--) { sum += *rp++ * *rp2++; } cs1 += *rp++ * *rp2++; /* right column sum */ } partsppmm[1][0] = cs0; partsppmm[1][1] = sum; partsppmm[1][2] = cs1; rp = row; rp2 = rp+nxs-1; partsppmm[2][0] = *rp++ * *rp2++; i = nx-3; sum = 0.0; while (i--) { sum += *rp++ * *rp2++; } partsppmm[2][1] = sum; partsppmm[2][2] = *rp++ * *rp2++; /* now each sum of cross terms is a sum of 4 parts */ for (j=0;j<2;j++) { for (i=0;i<2;i++) { xppmm[j][i] = 0.0; for (jj=0;jj<2;jj++) { for (ii=0;ii<2;ii++) { xppmm[j][i] += partsppmm[j+jj][i+ii]; }}}} /* next the cases with values from both images */ cmm=c0m=cpm=cm0=c00=cp0=cmp=c0p=cpp=sumxx = 0.0; nym2 = ny-2; j = nym2; nxm2 = nx-2; row = r; rowq = x+1+nxs; while (j--) { rp = row; rp2 = rp + 1; rp3 = rp +2; qp = rowq; row += nxs; rowq += nxs; i = nxm2; t2 = *rp++; t3 = *rp++; while (i--) { t0 = *qp++; t1 = t2; t2 = t3; t3 = *rp++; sumxx += t0 * t0; cmm += t0 * t1; c0m += t0 * t2; cpm += t0 * t3; } } row = r+nxs; rowq = x+1+nxs; j = nym2; while (j--) { rp = row; rp2 = rp + 1; rp3 = rp +2; qp = rowq; row += nxs; rowq += nxs; i = nxm2; t2 = *rp++; t3 = *rp++; while (i--) { t0 = *qp++; t1 = t2; t2 = t3; t3 = *rp++; cm0 += t0 * t1; c00 += t0 * t2; cp0 += t0 * t3; } } row = r+nxs+nxs; rowq = x+1+nxs; j = nym2; while (j--) { rp = row; rp2 = rp + 1; rp3 = rp +2; qp = rowq; row += nxs; rowq += nxs; i = nxm2; t2 = *rp++; t3 = *rp++; while (i--) { t0 = *qp++; t1 = t2; t2 = t3; t3 = *rp++; cmp += t0 * t1; c0p += t0 * t2; cpp += t0 * t3; } } /* have all the terms, now collect ones needed for quad 1 */ a1 = xx[1][1]; a4 = xx[0][0]; a2 = xx[1][0]; a3 = xx[0][1]; a5 = 2.0*xdx[1][0]; a10 = 2.0*xdx[0][0]; a9 = 2.0*xdy[0][0]; a6 = 2.0*xdy[0][1]; a7 = 2.0*xmmpp[0][0]; a8 = 2.0*xppmm[0][0]; a15 = sumxx; a11 = -2.*c00; a12 = -2.*cm0; a13 = -2.*c0m; a14 = -2.*cmm; /* printf("QUAD 1\n"); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10,a11 = %14.6e %14.6e %14.6e %14.6e\n",a8,a9,a10,a11); printf("a12,a13,a14,a15 = %14.6e %14.6e %14.6e %14.6e\n",a12,a13,a14,a15); */ mflag = 0; /* counts the number of mins */ qbest = 0.0; /* the best */ outside = 0.0; xoff = yoff = 0.0; getsxsy(); /* check if we got a min, only if somebody not zero */ if (subdx != 0 || subdy != 0) { xoff = - subdx; yoff = -subdy; qbest = mert(subdx, subdy); mflag += 1; if (subdx<0) outside = -subdx; else if (subdx>1.0) outside = subdx-1.0; if (subdy<0) outside = MAX(-subdy, outside); else if (subdy>1.0) outside = MAX(outside, subdy-1.0); } a4 = xx[0][2]; a2 = xx[1][2]; a5 = 2.0*xdx[1][1]; a10 = 2.0*xdx[0][1]; a9 = 2.0*xdy[0][2]; a7 = 2.0*xppmm[0][1]; a8 = 2.0*xmmpp[0][1]; a12 = -2.*cp0; a14 = -2.*cpm; /* printf("QUAD 2\n"); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10,a11 = %14.6e %14.6e %14.6e %14.6e\n",a8,a9,a10,a11); printf("a12,a13,a14,a15 = %14.6e %14.6e %14.6e %14.6e\n",a12,a13,a14,a15); */ getsxsy(); /* check if we got a min, only if somebody not zero */ if (subdx != 0 || subdy != 0) { qcur = mert(subdx, subdy); /* find our outside status */ qd = 0.0; if (subdx<0) qd = -subdx; else if (subdx>1.0) qd = subdx-1.0; if (subdy<0) qd = MAX(-subdy, qd); else if (subdy>1.0) qd = MAX(qd, subdy-1.0); if ( mflag == 0 ) { xoff = subdx; yoff = -subdy; qbest = qcur; outside = qd; } else { /* we have competition */ if ( outside != 0.0) { /* we are a contender, we decide on the outside issue if the previous was outside without reference to merit */ if (qd < outside) { outside = qd; xoff = subdx; yoff = -subdy; qbest = qcur;} } else { /* the competition has an outside of zip, if we also have zip then we compare on merit, if we have an outside problem here then the previous guy gets the honors */ if (qd == 0.0 && (qbest > qcur)) {xoff = subdx; yoff = -subdy; qbest = qcur;}} } mflag += 1; } a4 = xx[2][2]; a3 = xx[2][1]; a10 = 2.0*xdx[2][1]; a9 = 2.0*xdy[1][2]; a6 = 2.0*xdy[1][1]; a7 = 2.0*xmmpp[1][1]; a8 = 2.0*xppmm[1][1]; a13 = -2.*c0p; a14 = -2.*cpp; /* printf("QUAD 4\n"); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10,a11 = %14.6e %14.6e %14.6e %14.6e\n",a8,a9,a10,a11); printf("a12,a13,a14,a15 = %14.6e %14.6e %14.6e %14.6e\n",a12,a13,a14,a15); */ getsxsy(); /* check if we got a min, only if somebody not zero */ if (subdx != 0 || subdy != 0) { qcur = mert(subdx, subdy); /* find our outside status */ qd = 0.0; if (subdx<0) qd = -subdx; else if (subdx>1.0) qd = subdx-1.0; if (subdy<0) qd = MAX(-subdy, qd); else if (subdy>1.0) qd = MAX(subdy-1.0, qd); if ( mflag == 0 ) { xoff = subdx; yoff = subdy; qbest = qcur; outside = qd; } else { /* we have competition */ if ( outside != 0.0) { /* we are a contender, we decide on the outside issue if the previous was outside without reference to merit */ if (qd < outside) { outside = qd; xoff = subdx; yoff = subdy; qbest = qcur;} } else { /* the competition has an outside of zip, if we also have zip then we compare on merit, if we have an outside problem here then the previous guy gets the honors */ if (qd == 0.0 && (qbest > qcur)) {xoff = subdx; yoff = subdy; qbest = qcur;}} } mflag += 1; } a4 = xx[2][0]; a2 = xx[1][0]; a5 = 2.0*xdx[1][0]; a10 = 2.0*xdx[2][0]; a9 = 2.0*xdy[1][0]; a7 = 2.0*xppmm[1][0]; a8 = 2.0*xmmpp[1][0]; a12 = -2.*cm0; a14 = -2.*cmp; /* printf("QUAD 3\n"); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10,a11 = %14.6e %14.6e %14.6e %14.6e\n",a8,a9,a10,a11); printf("a12,a13,a14,a15 = %14.6e %14.6e %14.6e %14.6e\n",a12,a13,a14,a15); */ getsxsy(); /* check if we got a min, only if somebody not zero */ if (subdx != 0 || subdy != 0) { qcur = mert(subdx, subdy); /* find our outside status */ qd = 0.0; if (subdx<0) qd = -subdx; else if (subdx>1.0) qd = subdx-1.0; if (subdy<0) qd = MAX(-subdy, qd); else if (subdy>1.0) qd = MAX(qd, subdy-1.0); if ( mflag == 0 ) { xoff = -subdx; yoff = subdy; qbest = qcur; outside = qd; } else { /* we have competition */ if ( outside != 0.0) { /* we are a contender, we decide on the outside issue if the previous was outside without reference to merit */ if (qd < outside) { outside = qd; xoff = -subdx; yoff = subdy; qbest = qcur;} } else { /* the competition has an outside of zip, if we also have zip then we compare on merit, if we have an outside problem here then the previous guy gets the honors */ if (qd == 0.0 && (qbest > qcur)) {xoff = -subdx; yoff = subdy; qbest = qcur;}} } mflag += 1; } /* mflag contains the number of mins in case we want to document that or use the no min condition to search further */ //if (mflag == 0) printf("no min found\n"); if (mflag == 0) { xoff = yoff = 0.0; } /* zero these for no min case */ return; } /*------------------------------------------------------------------------- */ double sxvaluec(sy) double sy; { double c1, c2, xq; /* 3/4/97 changed a7+a8 to just a7 */ c1 = a5 + sy*a7 + sy*sy*a10; c2 = a2 + sy*sy*a4 +sy*a9; if (c2 == 0.0) { /* bump badmatch instead of printing */ /* printf("sxvaluec singular, c1, c2 = %g %g\n", c1,c2); */ badmatch++; return 0.0;} if ( isnan(c1) == 1 || isnan(c2) ==1 ) { printf("sxvaluec problem:\n"); printf("sy = %14.6e %14.6e\n", sy); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10 = %14.6e %14.6e %14.6e\n",a8,a9,a10); return 0.0; } xq = -.5*c1/c2; return xq; } /*------------------------------------------------------------------------- */ double syvaluec(sx) double sx; { double c1, c2, xq; /* 3/4/97 changed a7+a8 to just a7 */ c1 = a6 + sx*a7 + sx*sx*a9; c2 = a3 + sx*sx*a4 +sx*a10; if (c2 == 0.0) { /* bump badmatch instead of printing */ /* printf("syvalue singular, c1, c2 = %g %g\n",c1,c2); */ badmatch++; return 0.0;} if ( isnan(c1) == 1 || isnan(c2) ==1 ) { printf("syvaluec problem:\n"); printf("sx = %14.6e %14.6e\n", sx); printf("a1,a2,a3,a4 = %14.6e %14.6e %14.6e %14.6e\n",a1,a2,a3,a4); printf("a5,a6,a7 = %14.6e %14.6e %14.6e \n",a5,a6,a7); printf("a8,a9,a10 = %14.6e %14.6e %14.6e\n",a8,a9,a10); return 0.0; } xq = -.5*c1/c2; return xq; } /*------------------------------------------------------------------------- */ double subshiftc(xa, xb, nx, ny) /* we assume that the arrays have already been converted to F*8 */ double *xa, *xb; int nx, ny; { int nxs, nys; double t2, t1, t4, t3, d1, d2, d3, d4, sxz, syz; double x0, x1, x2, x3, y0, y1, y2, y3; int i, j, ii, jj, n, stride; double *xpa1, *xpa2, *xpb1, *xpb2; nxs = nx; nys = ny; stride = nxs - nx; a1=a2=a3=a4=a5=a6=a7=a8=a9=a10=0.0; xpa1 = xa; xpa2 = xa + nxs; xpb1 = xb; xpb2 = xb + nxs; j = ny - 1; while (j--) { x0 = *xpa1++; y0 = *xpb1++; x2 = *xpa2++; y2 = *xpb2++; i = nx -1; while (i--) { y3 = *xpb2++; x1 = *xpa1++; y1 = *xpb1++; x3 = *xpa2++; d1 = x0 - y3; d2 = x2 - y1; d3 = x1 - y2; d4 = x3 - y0; t1 = d1 + d2 + d3 + d4; t2 = d1 + d2 - d3 - d4; t3 = d1 - d2 + d3 - d4; t4 = d1 - d2 - d3 + d4; x0 = x1; x2 = x3; y0 = y1; y2 = y3; a1 += t1*t1; a2 += t2*t2; a3 += t3*t3; a4 += t4*t4; a5 += t1*t2; a6 += t1*t3; a7 += t1*t4; a8 += t2*t3; a9 += t2*t4; a10 += t3*t4; } xpa1 += stride; xpa2 += stride; xpb1 += stride; xpb2 += stride; } a1 = .0625*a1; a2 = .25 * a2; a3 = .25 * a3; a5 = .25 * a5; a6 = .25 * a6; a7 = .5 * a7; a8 = .5 * a8; /* 3/4/97 changed a7 to be a7+a8 */ a7 = a7 + a8; /* only 1 quad here so do the max find in line */ xoff = yoff = 0.0; n = 11; syz = 0.0; /* zero seed */ while (n--) { /* get new sxz */ sxz = sxvaluec(syz); if ( sxz < -0.5 || sxz > 0.5) { badmatch++; xoff = yoff = 0.0; return; } /* and get a new syz */ syz = syvaluec(sxz); if ( syz < -0.5 || syz > 0.5) { badmatch++; xoff = yoff = 0.0; return; } /* printf(" sxz, syz = %g %g\n", sxz,syz); */ } xoff = 2.*sxz; yoff = 2.*syz; return mertc(sxz, syz); } /*------------------------------------------------------------------------- */ double subshiftc_apod(xa, xb, gg, nx, ny) /* this version includes an apodizing function already prepared in gg which must be dimensioned (nx-1) by (ny-1) */ /* we assume that the arrays have already been converted to F*8 */ double *xa, *xb, *gg; int nx, ny; { int nxs, nys; double t2, t1, t4, t3, d1, d2, d3, d4, sxz, syz; double x0, x1, x2, x3, y0, y1, y2, y3, gapod; int i, j, ii, jj, n, stride; double *xpa1, *xpa2, *xpb1, *xpb2, xq; nxs = nx; nys = ny; stride = nxs - nx; a1=a2=a3=a4=a5=a6=a7=a8=a9=a10=0.0; xpa1 = xa; xpa2 = xa + nxs; xpb1 = xb; xpb2 = xb + nxs; j = ny - 1; while (j--) { x0 = *xpa1++; y0 = *xpb1++; x2 = *xpa2++; y2 = *xpb2++; i = nx -1; while (i--) { y3 = *xpb2++; x1 = *xpa1++; y1 = *xpb1++; x3 = *xpa2++; gapod = *gg++; d1 = x0 - y3; d2 = x2 - y1; d3 = x1 - y2; d4 = x3 - y0; t1 = d1 + d2 + d3 + d4; t2 = d1 + d2 - d3 - d4; t3 = d1 - d2 + d3 - d4; t4 = d1 - d2 - d3 + d4; x0 = x1; x2 = x3; y0 = y1; y2 = y3; xq = t1*gapod; a1 += xq*t1; a5 += xq*t2; a6 += xq*t3; a7 += xq*t4; xq = t2*gapod; a2 += xq*t2; a8 += xq*t3; a9 += xq*t4; xq = t3*gapod; a3 += xq*t3; a10 += xq*t4; a4 += gapod*t4*t4; } xpa1 += stride; xpa2 += stride; xpb1 += stride; xpb2 += stride; } a1 = .0625*a1; a2 = .25 * a2; a3 = .25 * a3; a5 = .25 * a5; a6 = .25 * a6; a7 = .5 * a7; a8 = .5 * a8; /* 3/4/97 changed a7 to be a7+a8 */ a7 = a7 + a8; /* only 1 quad here so do the max find in line */ xoff = yoff = 0.0; n = 11; syz = 0.0; /* zero seed */ while (n--) { /* get new sxz */ sxz = sxvaluec(syz); if ( sxz < -0.5 || sxz > 0.5) { xoff = yoff = 0.0; return; } /* and get a new syz */ syz = syvaluec(sxz); if ( syz < -0.5 || syz > 0.5) { xoff = yoff = 0.0; return; } /* printf(" sxz, syz = %g %g\n", sxz,syz); */ } xoff = 2.*sxz; yoff = 2.*syz; return mertc(sxz, syz); } /*------------------------------------------------------------------------- */