include include include include include include define RSV_LABLEN 60 # RDIFFERX -- program to produce histograms of reseau mark shifts in # various ways # # D. Giaretta, 1-Aug-1987 Original SPP version procedure t_rdifferx() char input1[SZ_FNAME] # entry file - without whitespace char entry[SZ_R_PATTERN] # first entry/pattern char input2[SZ_FNAME] # reference file - without whitespace char refentry[SZ_R_PATTERN] # reference entry char output[SZ_FNAME] # output reseau file char option[SZ_LINE] # plotoption int nbins # number of bins real z1, z2 # hist max and min - if INDEF then # use data max and min bool listout # just print data instead of plotting char device[SZ_FNAME] # graphics device real diffmean # mean radial diff real diffmax # max radial diff real diffrms # rms radial diff int ngood # number of good points int nrows # number of rows in grid int ncols # number of cols in grid #-- pointer rp # pointer to input1 pointer rentp # pointer to entry data pointer rpref # pointers for ref. file pointer refpt # pointer to ref data pointer rpout # pointer to output file pointer routpt # pointer to output data pointer tempp char reply[1], stat, getchar() # for delays before frame bool reffile # was separate ref file given? real xlims[2, 3] # limits int i bool clgetb(), first, optall int clgeti() pointer dp, hp # structures for diff data and # histograms pointer gp, gopen() # graphics stream real clgetr() int kwindex, clgwrd() pointer rs_template(), rlist, rs_open() int rs_readr(), rs_list(), nowhite(), junk bool rs_check(), streq() errchk rs_open, rs_template, rs_readr, rs_list errchk gopen begin rpref = NULL # get input1 and check for blanks call clgstr( "input1", input1, SZ_FNAME) if ( nowhite( input1, input1, SZ_FNAME) == 0 ) call error( 0, "input file must be given") call clgstr( "entry", entry, SZ_R_PATTERN) call clgstr( "input2", input2, SZ_FNAME) if ( nowhite( input2, input2, SZ_FNAME) == 0 ) reffile = false else reffile = true # get entry defs call clgstr( "refentry", refentry, SZ_R_PATTERN) call clgstr( "output", output, SZ_FNAME) junk = nowhite( output, output, SZ_FNAME) # open input1 and output - may be the same if ( streq( input1, output) ) { rp = rs_open( input1, READ_WRITE, NULL) rpout = rp } else { rp = rs_open( input1, READ_ONLY, NULL) if ( streq( input2, output ) ) { rpout = rs_open( output, READ_WRITE, NULL) rpref = rpout } else { rpout = rs_open( output, NEW_COPY, rp) } } # get options - check for all vectors on one plot kwindex = clgwrd( "option", option, SZ_LINE, ",all,sequential,ref") optall = false if ( kwindex == 1 ) optall = true else if ( kwindex == 2 && reffile ) # but we ignore input2 if sequential pairs reffile = false # input2 if required if (reffile && rpref == NULL ){ rpref = rs_open( input2, READ_ONLY, 0) # check if grid sizes the same if (!rs_check(rpref, RES_NROWS(rp), RES_NCOLS(rp)) ) call error (0, "input1 and input2 grid sizes must match") } else { # if no input2, set ref. pointer to input1 rpref = rp call strcpy( input1, input2, SZ_FNAME) } nbins = clgeti( "nbins") listout = clgetb( "listout") z1 = clgetr( "z1") z2 = clgetr( "z2") for (i=1; i<=3; i=i+1){ xlims[1, i] = z1 xlims[2, i] = z2 } call clgstr("device", device, SZ_FNAME] # open the device, if required if ( !listout) { gp = gopen(device, NEW_FILE, STDGRAPH) # set the graphics stuff call gsetr( gp, G_TXSIZE, 0.5) call gseti( gp, G_NTITLELINES, 0) call gsetr( gp, G_ASPECT, 1.0) } # allocate space for reference call rs_alloc( rp, TY_REAL, rentp ) call rs_alloc( rpref, TY_REAL, refpt ) call rs_alloc( rp, TY_REAL, routpt) # allocate space for histogram call salloc( dp, RES_EXTENT(rp), TY_REAL) call salloc( hp, nbins+2, TY_REAL) # read reference entry rlist = rs_template( rpref, refentry) if ( rs_readr( rpref, refpt, EOS, rs_list(rlist) ) == RES_F_NONEXTENTRY) call error( 0, "could not read reference entry") # initialise entry loop rlist = rs_template( rp, entry) first = true # plot until we reach the end of list while (rs_readr( rp, rentp, EOS, rs_list(rlist)) != RES_F_NONEXTENTRY) { # wait for keystroke before starting next screen, if we are # going to plot more frames, but omit first time if (!first && !optall && !listout){ call printf( "press RETURN for next plot:") call flush(STDOUT) stat = getchar(reply) call gframe(gp) } call xr_hist( gp, input1, input2, first, optall, entry, nbins, z1, z2, listout, option, rp, rentp, rpref, refpt, Memr[RES_XPT(rentp)], Memr[RES_YPT(rentp)], Memb[RES_FPT(rentp)], Memr[RES_XPT(refpt)], Memr[RES_YPT(refpt)], Memb[RES_FPT(refpt)], Memr[dp], Memr[hp], xlims, routpt, Memr[RES_XPT(routpt)], Memr[RES_YPT(routpt)], Memb[RES_FPT(routpt)] ) # if we are plotting sequential pairs we must swop the # pointers if (option[1] == 's' || option[1] == 'S'){ tempp = refpt refpt = rentp rentp = tempp } if ( rpout != NULL) call rs_write( rpout, routpt, 0) first = false } # calculate the output parameters call rsd_para( routpt, Memr[RES_XPT(routpt)], Memr[RES_YPT(routpt)], Memb[RES_FPT(routpt)], diffmean, diffrms, diffmax, nrows, ncols, ngood) call clputr( "diffmean", diffmean) call clputr( "diffmax", diffmax) call clputr( "diffrms", diffrms) call clputi( "nrows", nrows ) call clputi( "ncols", ncols ) call clputi( "ngood", ngood ) # remember we have swopped refpt and rentp call clpstr( "refent", RES_ENT_ENTRY(rentp), SZ_R_ENTRY) call clpstr( "resent", RES_ENT_ENTRY(refpt), SZ_R_ENTRY) # close the graphics stream call gclose(gp) # close the files call rs_close(rp) call rs_close(rpref) call rs_close(rpout) end # XR_HIST -- plot reseau shift histograms procedure xr_hist( gp, input1, input2, first, optall, entry, nbins, z1, z2, listout, option, rp, rentp, rpref, refpt, x, y, f, xref, yref, fref, data, hist, xlims, routpt, xout, yout, fout ) pointer gp # i: graphics stream pointer char input1[SZ_FNAME] # i: name of input1 char input2[SZ_FNAME] # i: name of input2 bool first, optall # i: flags char entry[SZ_R_PATTERN] # i: entry names int nbins # i: number of bins real z1, z2 # i: limits for plots bool listout # i: list only? char option[SZ_LINE] # i: plotting option pointer rp, rpref # i: reference reseau pointers pointer rentp, refpt # i: reseau pointers real x[ARB], xref[ARB] # i: x coords real y[ARB], yref[ARB] # i: y coords bool f[ARB], fref[ARB] # i: null res. ? flags # true if bad point real data[ARB] # o: for storing diffs real hist[ARB] # o: to store histogram data real xlims[2, 3] # i/o: x axis limits pointer routpt # i: pointer for output entry real xout[ARB] # o: output x data real yout[ARB] # o: output y data real fout[ARB] # o: output flag data #-- char date[SZ_TIME], datetime[SZ_R_TIME] long ldate int nrows, ncols char msg1[RSV_LABLEN], msg2[RSV_LABLEN] char msg3[RSV_LABLEN], msg4[RSV_LABLEN] char msg5[RSV_LABLEN] long clktime() begin ncols = RES_ENT_NCOLS(rentp) nrows = RES_ENT_NROWS(rentp) # calculate the histograms # plot if required if ( !listout) { # if we as just drawing first plot, or if we are not plotting # everything on one plot, initialise the window if (!optall || first){ # set up headers for printing call sprintf( msg1, RSV_LABLEN, "diff file : %s \n") call pargstr( input1 ) if (!optall){ call rs_cnvtime( rentp, datetime, SZ_R_TIME) call sprintf( msg2, RSV_LABLEN, "run entry : %s , dated %s \n") call pargstr( RES_ENT_ENTRY(rentp) ) call pargstr( datetime ) } else { call sprintf( msg2, RSV_LABLEN, "run entries: %s \n") call pargstr( entry ) } call sprintf( msg3, RSV_LABLEN, "ref file : %s \n") call pargstr( input2 ) call rs_cnvtime( refpt, datetime, SZ_R_TIME) call sprintf( msg4, RSV_LABLEN, "ref entry : %s , dated %s \n") call pargstr( RES_ENT_ENTRY(refpt) ) call pargstr( datetime ) ldate = clktime(0) call cnvtime(ldate, date, SZ_TIME] call sprintf( msg5, RSV_LABLEN, "plotted : %s\n") call pargstr( date ) # write labels call gseti( gp, G_WCS, 0) call gtext( gp, 0.1, 0.97, msg1, EOS) call gtext( gp, 0.1, 0.94, msg2, EOS) call gtext( gp, 0.1, 0.91, msg3, EOS) call gtext( gp, 0.1, 0.88, msg4, EOS) call gtext( gp, 0.1, 0.85, msg5, EOS) } } # plot or print differences call xrd_hist( gp, listout, nbins, z1, z2, rp, rentp, rpref, refpt, x, y, f, xref, yref, fref, data, hist, xlims, xout, yout, fout ) # copy header entry info if required RES_ENT_NROWS(routpt) = RES_ENT_NROWS(rentp) RES_ENT_NCOLS(routpt) = RES_ENT_NCOLS(rentp) RES_ENT_SFIELD1(routpt) = RES_ENT_SFIELD1(rentp) RES_ENT_SFIELD2(routpt) = RES_ENT_SFIELD2(rentp) RES_ENT_SFIELD3(routpt) = RES_ENT_SFIELD3(rentp) RES_ENT_SFIELD4(routpt) = RES_ENT_SFIELD4(rentp) RES_ENT_DATE(rentp) = RES_ENT_DATE(routpt) call strcpy(RES_ENT_ENTRY(rentp), RES_ENT_ENTRY(routpt), SZ_R_ENTRY) call strcpy(RES_ENT_TRACKING(rentp), RES_ENT_TRACKING(routpt), SZ_R_TRACKING) if (!listout) call gflush(gp) end # XRD_HIST -- plot histogram procedure xrd_hist( gp, listout, nbins, z1, z2, rp, rentp, rpref, refpt, x, y, f, xref, yref, fref, data, hist, xlims, xout, yout, fout ) pointer gp # i: graphics pointer bool listout # i: list only? int nbins # i: no. of bins real z1, z2 # i/o: limits pointer rp, rentp, rpref, refpt # I; reseau pointers real x[ARB], y[ARB] # i: x, y coords real xref[ARB], yref[ARB] # i: x, y ref coords bool f[ARB], fref[ARB] # i: null flags real data[ARB], hist[ARB] # o: diff and hist data real xlims[2, 3] # i/o: hist limits real xout[ARB] # o: output x data real yout[ARB] # o: output y data bool fout[ARB] # o: output flag data #-- int npts int i, pos real maxdiff, mindiff, width errchk gseti, gsview begin npts = RES_ENT_NROWS(rentp)*RES_ENT_NCOLS(rentp) # calculate output for ( i=1; i<=npts; i=i+1) { if ( !f[i] && !fref[i] ) { xout[i] = x[i] - xref[i] yout[i] = y[i] - yref[i] fout[i] = false } else { xout[i] = INDEFR yout[i] = INDEFR fout[i] = true } } # radial differences maxdiff = 0.0 mindiff = 0.0 for (i=1; i<= npts ;i=i+1) if (!f[i] && !fref[i]) { data[i] = sqrt((x[i]-xref[i])**2 + (y[i]-yref[i])**2) maxdiff = max( data[i], maxdiff) mindiff = min( data[i], mindiff) } else data[i] = INDEF if (xlims[1, 1] == INDEF) xlims[1, 1] = mindiff if (xlims[2, 1] == INDEF) xlims[2, 1] = maxdiff width = max( (xlims[2, 1] - xlims[1, 1])/(nbins-1), 1.0e-7) for (i=1; i<=nbins; i=i+1) hist[i] = 0 for (i=1; i<=npts; i=i+1) { if (!IS_INDEFR (data[i])) { pos = int((data[i] - xlims[1, 1])/width) + 1 if (pos < 1) pos = 1 else if (pos >nbins) pos = nbins hist[pos] = hist[pos] + 1 } } if (!listout) { call gseti(gp, G_WCS, 1) call gsview( gp, 0.2, 1.0, 0.6, 0.85) call xr_histcrv( gp, nbins, npts, data, hist, xlims[1, 1], EOS, EOS, "Rad diffs") } else { for (i=1; i<=nbins ; i=i+1) { call printf( " %12d %6d \n") call pargi( int( hist[i] ) ) call pargi( i) } } # x direction differences maxdiff = -MAX_REAL mindiff = +MAX_REAL for (i=1; i<= npts ;i=i+1) if (!f[i] && !fref[i]) { data[i] = x[i]-xref[i] maxdiff = max(data[i], maxdiff) mindiff = min(data[i], mindiff) } else data[i] = INDEF if (xlims[1, 2] == INDEF) xlims[1, 2] = mindiff if (xlims[2, 2] == INDEF) xlims[2, 2] = maxdiff width = max( (xlims[2, 2] - xlims[1, 2])/(nbins-1), 1.0e-7) for (i=1; i<=nbins; i=i+1) hist[i] = 0 for (i=1; i<=npts; i=i+1) { if (!IS_INDEFR (data[i])) pos = int((data[i] - xlims[1, 2])/width) + 1 if (pos < 1) pos = 1 else if (pos >nbins) pos = nbins hist[pos] = hist[pos] + 1 } if (!listout) { call gseti(gp, G_WCS, 2) call gsview( gp, 0.2, 1.0, 0.35, 0.6) call xr_histcrv( gp, nbins, npts, data, hist, xlims[1, 2], EOS, EOS, "X diffs") } else { for (i=1; i<=nbins ; i=i+1) { call printf( " %12d %6d \n") call pargi( int( hist[i] ) ) call pargi( i) } } # y direction differences maxdiff = -MAX_REAL mindiff = +MAX_REAL for (i=1; i<= npts ;i=i+1) if (!f[i] && !fref[i]) { data[i] = y[i]-yref[i] maxdiff = max(data[i], maxdiff) mindiff = min(data[i], mindiff) } else data[i] = INDEF if (xlims[1, 3] == INDEF) xlims[1, 3] = mindiff if (xlims[2, 3] == INDEF) xlims[2, 3] = maxdiff width = max( (xlims[2, 3] - xlims[1, 3])/(nbins-1), 1.0e-7) for (i=1; i<=nbins; i=i+1) hist[i] = 0 for (i=1; i<=npts; i=i+1) { if (!IS_INDEFR (data[i])) pos = int((data[i] - xlims[1, 3])/width) + 1 if (pos < 1) pos = 1 else if (pos >nbins) pos = nbins hist[pos] = hist[pos] + 1 } if (!listout) { call gseti(gp, G_WCS, 3) call gsview( gp, 0.2, 1.0, 0.1, 0.35) call xr_histcrv( gp, nbins, npts, data, hist, xlims[1, 3], EOS, EOS, "Y diffs") } else { for (i=1; i<=nbins ; i=i+1) { call printf( " %12d %6d \n") call pargi( int( hist[i] ) ) call pargi( i) } } end # XR_HISTCRV -- plot data as histogram procedure xr_histcrv( gp, nbins, npts, data, hist, xlims, title, xlab, ylab) pointer gp # i: Graphics descriptor int nbins # i: no. of bins int npts # i: number of data points real data[ARB] # i: Vector to be plotted (Y values) real hist[ARB] real xlims[2] char title[ARB] char xlab[ARB] char ylab[ARB] #-- errchk gswind, glabax, xhstcrv begin call gswind( gp, xlims[1], xlims[2], 0.0, float(npts)) call glabax (gp, title, xlab, ylab) call xhstcrv (gp, hist, nbins, xlims[1], xlims[2], GF_HOLLOW) end # XHSTCRV -- Draws a histogram style curve (bar graph). procedure xhstcrv (gp, v, npts, xmin, xmax, style) pointer gp # i: Graphics descriptor real v[ARB] # i: Vector to be plotted (Y values) int npts # i: Number of data values real xmin, xmax # i: Range of vector in X int style # i: Fill area style int bin real x real xbin[4], ybin[4] # Fill polyline real dx, bw real left, right, bottom, top # Data window errchk ggwind, gfill begin dx = (xmax - xmin) / (npts - 1) call ggwind (gp, left, right, bottom, top) # Bin width is in world coordinates bw = -dx/2.0 do bin = 1, npts { if (!IS_INDEF(v[bin])) { x = xmin + dx*real (bin-1) xbin[1] = x - bw ybin[1] = bottom xbin[2] = xbin[1] ybin[2] = v[bin] xbin[3] = x + bw ybin[3] = v[bin] xbin[4] = xbin[3] ybin[4] = bottom call gfill (gp, xbin, ybin, 4, style) } } end # RSD_PARA -- calculate the gross differences for the last pair procedure rsd_para( routpt, x, y, f, diffmean, diffrms, diffmax, nrows, ncols, ngood) pointer routpt # i: pointer for entry real x[ARB] # i: x data real y[ARB] # i: y data bool f[ARB] # i: flag data real diffmean # o: mean difference real diffrms # o: rms difference real diffmax # o: max radial difference int nrows # o: number of rows int ncols # o: number of cols int ngood # o: number of good points #-- int i double sum2, sum, temp, maxtemp begin sum = 0.0 sum2 = 0.0 ngood = 0 nrows = RES_ENT_NROWS(routpt) ncols = RES_ENT_NCOLS(routpt) maxtemp = 0.0 for ( i=1; i<= nrows*ncols; i=i+1 ) { if ( !f[i] ) { ngood = ngood + 1 temp = x[i]**2 + y[i]**2 sum = sum + sqrt(temp) sum2 = sum2 + temp maxtemp = max( maxtemp, temp) } } if ( ngood == 0 ) { diffmean = INDEFR diffrms = INDEFR diffmax = INDEFR } else { diffmean = sum/float( ngood ) diffmax = sqrt(maxtemp) if ( ngood <= 1 ) diffrms = INDEFR else diffrms = sqrt ( sum2/float(ngood - 1) ) } end