diff options
Diffstat (limited to 'noao/imred/vtel/syndico.x')
| -rw-r--r-- | noao/imred/vtel/syndico.x | 416 |
1 files changed, 416 insertions, 0 deletions
diff --git a/noao/imred/vtel/syndico.x b/noao/imred/vtel/syndico.x new file mode 100644 index 00000000..64910679 --- /dev/null +++ b/noao/imred/vtel/syndico.x @@ -0,0 +1,416 @@ +include <mach.h> +include <imhdr.h> +include <imset.h> +include <gset.h> +include "syndico.h" +include "vt.h" + +# SYNDICO -- Make Dicomed prints of synoptic images. This program is tuned +# to make the images 18 centimeters in diameter. + +procedure t_syndico() + +char image[SZ_FNAME] # image to plot +char logofile[SZ_FNAME] # file containing logo +char device[SZ_FNAME] # plot device +int sbthresh # squibby brightness threshold +bool verbose # verbose flag +bool plotlogo # plotlogo flag +bool forcetype # force image type flag +bool magnetic # image type = magnetic flag + +int obsdate, wavelength, obstime +int i, j, month, day, year, hour, minute, second, stat, bufptr +real delta_gblock, x, y +real excen, eycen, exsmd, eysmd, rguess +real b0, l0 +real mapy1, mapy2, radius, scale, diskfrac +char ltext[SZ_LINE] +char system_id[SZ_LINE] + +short grey[16] +pointer gp, sp, im, lf +pointer subrasp, subras1, buff +int trnsfrm[513] +int lkup10830[1091] +int gs10830[16] +real xstart, xend, ystart, yend, yinc +real xcenerr, ycenerr, ndc_xcerr, ndc_ycerr +real temp_xcenter, temp_ycenter + +pointer immap(), gopen(), imgl2s() +int imgeti(), clgeti(), open(), read() +real imgetr() +bool clgetb(), imaccf() +include "trnsfrm.inc" +errchk gopen, immap, sysid, imgs2s, imgl2s + +# Grey scale points for 10830. +data (gs10830[i], i = 1, 6) /-1000,-700,-500,-400,-300,-250/ +data (gs10830[i], i = 7, 10) /-200,-150,-100,-50/ +data (gs10830[i], i = 11, 16) /0,10,20,40,60,90/ + +begin + call smark (sp) + call salloc (subrasp, DIM_VTFD, TY_SHORT) + call salloc (subras1, 185*185, TY_SHORT) + call salloc (buff, 185, TY_CHAR) + + # Get parameters from the cl. + call clgstr ("image", image, SZ_FNAME) + call clgstr ("logofile", logofile, SZ_FNAME) + call clgstr ("device", device, SZ_FNAME) + sbthresh = clgeti ("sbthresh") + plotlogo = clgetb ("plotlogo") + verbose = clgetb ("verbose") + forcetype = clgetb ("forcetype") + magnetic = clgetb ("magnetic") + + # Open the input image, open the logo image if requested. + im = immap (image, READ_ONLY, 0) + if (plotlogo) + iferr { + lf = open (logofile, READ_ONLY, TEXT_FILE) + } then { + call eprintf ("Error opening the logo file, logo not made.\n") + plotlogo = false + } + + # Get/calculate some of the housekeeping data. + if (imaccf (im, "obs_date")) { + obsdate = imgeti (im, "obs_date") + obstime = imgeti (im, "obs_time") + month = obsdate/10000 + day = obsdate/100 - 100 * (obsdate/10000) + year = obsdate - 100 * (obsdate/100) + hour = int(obstime/3600) + minute = int((obstime - hour * 3600)/60) + second = obstime - hour * 3600 - minute * 60 + } else { + # Use cl query parameters to get these values. + call eprintf ("Date and Time not found in image header.\n") + call eprintf ("Please enter them below.\n") + month = clgeti ("month") + day = clgeti ("day") + year = clgeti ("year") + hour = clgeti ("hour") + minute = clgeti ("minute") + second = clgeti ("second") + } + + # Get the solar image center and radius from the image header, + # get the solar image radius from the ephemeris routine. If + # the two radii are similar, use the former one, if they are + # %10 percent or more different, use the ephemeris radius and + # assume the center is at (1024,1024). + + # Get ellipse parameters from image header. + # If they are not there, warn the user that we are using ephemeris + # values. + if (imaccf (im, "E_XCEN")) { + excen = imgetr (im, "E_XCEN") + eycen = imgetr (im, "E_YCEN") + exsmd = imgetr (im, "E_XSMD") + eysmd = imgetr (im, "E_YSMD") + + # Get rguess from ephem. + iferr (call ephem (month, day, year, hour, minute, second, rguess, + b0, l0, false)) + call eprintf ("Error getting ephemeris data.\n") + + radius = (exsmd + eysmd) / 2.0 + if (abs(abs(radius-rguess)/rguess - 1.0) > 0.1) { + radius = rguess + excen = 1024.0 + eycen = 1024.0 + } + + } else { + call eprintf ("No ellipse parameters in image header.\n Using") + call eprintf (" ephemeris value for radius and setting center to") + call eprintf (" 1024, 1024\n") + + # Get rguess from ephem. + iferr (call ephem (month, day, year, hour, minute, second, rguess, + b0, l0, false)) + call eprintf ("Error getting ephemeris data.\n") + + radius = rguess + excen = 1024.0 + eycen = 1024.0 + } + + # Error in center. (units of pixels) + xcenerr = excen - 1024.0 + ycenerr = eycen - 1024.0 + + # Transform error to NDC. + ndc_xcerr = xcenerr * (1.0/4096.0) + ndc_ycerr = ycenerr * (1.0/4096.0) + + # Next, knowing that the image diameter must be 18 centimeters, + # calculate the scaling factor we must use to expand the image. + # DICO_18CM is a MAGIC number = 18 centimeters on dicomed prints + # given the way the NOAO photo lab currently enlarges the images. + scale = DICO_18CM / real(radius*2) + + # Open the output file. + gp = gopen (device, NEW_FILE, STDGRAPH) + + # Put feducial(sp?) marks on plot. + diskfrac = radius/1024.0 + temp_xcenter = DICO_XCENTER-ndc_xcerr + temp_ycenter = DICO_YCENTER-ndc_ycerr + call gline (gp, temp_xcenter, temp_ycenter+diskfrac*.25*scale+.01, + temp_xcenter, temp_ycenter+diskfrac*.25*scale+.025) + call gline (gp, temp_xcenter, temp_ycenter-diskfrac*.25*scale-.01, + temp_xcenter, temp_ycenter-diskfrac*.25*scale-.025) + + # Draw a little compass on the plot. + call gline (gp, .25, DICO_YCENTER+.25+.01, + .25, DICO_YCENTER+.25+.035) + call gtext (gp, .25, DICO_YCENTER+.25+.037, + "N", "v=b;h=c;s=.50") + call gmark (gp, .2565, DICO_YCENTER+.25+.037, + GM_CIRCLE, .006, .006) + call gmark (gp, .2565, DICO_YCENTER+.25+.037, + GM_CIRCLE, .001, .001) + call gline (gp, .25, DICO_YCENTER+.25+.01, + .28, DICO_YCENTER+.25+.01) + call gtext (gp, .282, DICO_YCENTER+.25+.01, + "W", "v=c;h=l;s=.50") + call gmark (gp, .290, DICO_YCENTER+.25+.01-.006, + GM_CIRCLE, .006, .006) + call gmark (gp, .290, DICO_YCENTER+.25+.01-.006, + GM_CIRCLE, .001, .001) + + # Get the wavelength from the image header. If the user wants + # to force the wavelength, do so. (this is used if the header + # information about wavelength is wrong.) + wavelength = imgeti (im, "wv_lngth") + if (forcetype) + if (magnetic) + wavelength = 8688 + else + wavelength = 10830 + + # Write the grey scale labels onto the plot. + delta_gblock = (IMGTR_X - IMGBL_X)/16. + y = IMGBL_Y - .005 + do i = 1, 16 { + x = IMGBL_X + real(i-1) * delta_gblock + delta_gblock/2. + call sprintf (ltext, SZ_LINE, "%d") + if (wavelength == 8688) + call pargi ((i-1)*(int((512./15.)+0.5))-256) + else if (wavelength == 10830) + call pargi (gs10830(i)) + call gtext (gp, x, y, ltext, "v=t;h=c;s=.20") + } + + # Label on grey scale. + call sprintf (ltext, SZ_LINE, "%s") + if (wavelength == 8688) + call pargstr ("gauss") + else if (wavelength == 10830) + call pargstr ("relative line strength") + call gtext (gp, DICO_XCENTER, (IMGBL_Y-.024), ltext, "v=c;h=c;s=.5") + + # Put the title on. + call sprintf (ltext, SZ_LINE, "%s") + if (wavelength == 8688) + call pargstr ("8688 MAGNETOGRAM") + else if (wavelength == 10830) + call pargstr ("10830 SPECTROHELIOGRAM") + else + call pargstr (" ") + call gtext (gp, DICO_XCENTER, .135, ltext, "v=c;h=c;s=.7") + + # If we don't have a logo to plot, write the data origin on the plot. + if (!plotlogo) { + + call sprintf (ltext, SZ_LINE, "%s") + call pargstr ("National") + call gtext (gp, .24, .155, ltext, "v=c;h=c;s=.7") + call sprintf (ltext, SZ_LINE, "%s") + call pargstr ("Solar") + call gtext (gp, .24, .135, ltext, "v=c;h=c;s=.7") + call sprintf (ltext, SZ_LINE, "%s") + call pargstr ("Observatory") + call gtext (gp, .24, .115, ltext, "v=c;h=c;s=.7") + } + + # Put month/day/year on plot. + call sprintf (ltext, SZ_LINE, "%02d/%02d/%02d") + call pargi (month) + call pargi (day) + call pargi (year) + call gtext (gp, .70, .175, ltext, "v=c;h=l;s=.5") + + # Put the hour:minute:second on plot. + call sprintf (ltext, SZ_LINE, "%02d:%02d:%02d UT") + call pargi (hour) + call pargi (minute) + call pargi (second) + call gtext (gp, .70, .155, ltext, "v=c;h=l;s=.5") + + # Fill in the grey scale. + if (wavelength == 8688) { + do i = 1, 16 + grey[i] = (trnsfrm[(i-1)*(int((512./15.)+0.5))+1]) + call gpcell (gp, grey, 16, 1, IMGBL_X, IMGBL_Y, IMGTR_X, IMGTR_Y) + } else if (wavelength == 10830) { + do i = 1, 16 + grey[i] = (lkup10830[gs10830(i)+1001]) + call gpcell (gp, grey, 16, 1, IMGBL_X, IMGBL_Y, IMGTR_X, IMGTR_Y) + } + + # Prepare some constants for plotting. + xstart = temp_xcenter - .25 * scale + xend = temp_xcenter + .25 * scale + ystart = temp_ycenter - .25 * scale + yend = temp_ycenter + .5 * scale + mapy1 = ystart + mapy2 = ystart + yinc = (.5*scale)/real(DIM_VTFD) + + # Put the data on the plot. Line by line. + do i = 1, DIM_VTFD { + + if (verbose) { + call printf ("line = %d\n") + call pargi (i) + call flush (STDOUT) + } + + subrasp = imgl2s (im, i) + + # Call the limb trimmer and data divider. + call fixline (Mems[subrasp], DIM_VTFD, wavelength, sbthresh) + + # Update the top and bottom edges of this line. + mapy1 = mapy2 + mapy2 = mapy2 + yinc + + # Put the line on the output plot. + call gpcell (gp, Mems[subrasp], DIM_VTFD, 1, xstart, + mapy1, xend, mapy2) + + } # End of do loop on image lines. + + # Put the system identification on the plot. + call sysid (system_id, SZ_LINE) + call gtext (gp, DICO_XCENTER, .076, system_id, "h=c;s=0.45") + + # Put the NSO logo on the plot. + if (plotlogo) { + + # Read in the image. (the image is encoded in a text file) + do i = 1, 185 { + bufptr = 0 + while (bufptr < 185-79) { + stat = read (lf, Memc[buff+bufptr], 80) + bufptr = bufptr + 79 + } + stat = read (lf, Memc[buff+bufptr], 80) + do j = 1, 185 { + Mems[subras1+(i-1)*185+j-1] = + short((Memc[buff+j-1]-32.)*2.7027027) + } + } + + # Put it on the plot. + call gpcell (gp, Mems[subras1], 185, 185, .24, .13, .32, .21) + } + + # Close the graphics pointer, unmap images, free stack. + call gclose (gp) + call imunmap (im) + if (plotlogo) + call close (lf) + call sfree (sp) +end + + +# FIXLINE -- Clean up the line. Set the value of pixels off the limb to +# zero, remove the squibby brightness from each pixel, and apply a +# nonlinear lookup table to the greyscale mapping. + +procedure fixline (ln, xlength, wavelength, sbthresh) + +int xlength # length of line buffer +short ln[xlength] # line buffer +int wavelength # wavelength of the observation +int sbthresh # squibby brightness threshold + +int trnsfrm[513] +int lkup10830[1091] +bool found +int i, left, right +include "trnsfrm.inc" + +begin + # Look in from the left end till squibby brightness goes above the + # threshold, remember where this limbpoint is. + found = false + do i = 1, xlength { # Find left limbpoint. + if (and(int(ln[i]),17B) > sbthresh) { + found = true + left = i + break + } + } + + if (found) { + # Find the right limbpoint. + do i = xlength, 1, -1 { + if (and(int(ln[i]),17B) > sbthresh) { + right = i + break + } + } + + # Divide the image by 16, map the greyscale, and trim the limb. + do i = left+1, right-1 { + + # Remove squibby brightness. + ln[i] = ln[i]/16 + + if (wavelength == 8688) { + # Magnetogram, nonlinear greyscale. + # Make data fit in the table. + if (ln[i] < -256) + ln[i] = -256 + if (ln[i] > 256) + ln[i] = 256 + + # Look it up in the table. + ln[i] = trnsfrm[ln[i]+257] + } else if (wavelength == 10830) { + # 10830 spectroheliogram, nonlinear greyscale. + # Make data fit in the table. + if (ln[i] < -1000) + ln[i] = -1000 + if (ln[i] > 90) + ln[i] = 90 + # Look it up in the table. + ln[i] = lkup10830[ln[i]+1001] + } else { + # Unknown type, linear greyscale. + if (ln[i] < 1) + ln[i] = 1 + if (ln[i] > 255) + ln[i] = 255 + } + } + + # Set stuff outside the limb to zero. + do i = 1, left + ln[i] = 0 + do i = right, xlength + ln[i] = 0 + } else { + # This line is off the limb, set it to zero. + do i = 1, xlength + ln[i] = 0 + } +end |