;------------------------------------------------------------------ ;+ ; NAME: ; SIMPMAP ; PURPOSE: ; This procedure is used to convert a BBSO VMG image to ; "Peter Wilson" format for the polar fields project. ; This procedure is a port of the 6/29/95 IRAF imfort ; task. ; ; The procedure takes a "squared" image and projects it ; onto a square grid of latitude vs. longitude. ; ; The procedure determines the latitude simply by taking the ; absolute value of the latitude it finds. Therefore, ; equatorial images, with both positive and negative latitudes, ; will not be correctly displayed. The program also ceneters ; its longitude display to either side of center sun. ; ; The logfile ("Peter Wilson format" file) produced by this ; procedure contains the image data output as a formatted file. ; The image is 181 rows by 181 columns; in other words, ; a grid in which each pixel is one degree of latitude by ; one degree of longitude. ; CATEGORY: ; POLAR ; CALLING SEQUENCE: ; SIMPMAP,imname,logfile ; INPUTS: ; imname = BBSO VMG image file to map. This image must be "squared"; ; that is, it must have square pixels. ; The actual coordinates of Sun center (h and k in pixels) ; must be in the FITS header keywords CENTERX and CENTERY. ; The radius of the Sun in pixels must be in the keyword ; RADIUS. The pixel coordinates are based on the lower ; left corner of the picture being (0,0). ; KEYWORD PARAMETERS: ; OUTPUTS: ; logfile = text file containing the map information output by ; the procedure. Each line output contains the ; longitude (zero at center), latitude, the total ; gray level at the pixel (summing the contribution ; from each image pixel), the number of image pixels ; contributing, and the final map value (the total ; divided by the number contributing). ; imnamem = FITS image (181 x 181 pixels) version of the final ; map. The name of the file is imname concatenated ; with 'm'. ; cimname = FITS image (181 x 181 pixels) containing the ; number of pixels contributing to the map at each point. ; The name of the file is 'c' concatenated with imname. ; COMMON BLOCKS: ; NOTES: ; The first 11 and last 8 columns in the VMGs are bad. ; The procedure does not use them and subtracts 11 from the ; h value before mapping. Uses the procedure sunlong to ; find P, B, and the Carrington longitude of the central meridian. ; MODIFICATION HISTORY: ; J. Varsik, 19 Oct. 1995 --- Port from IRAF task ;- ;------------------------------------------------------------------- pro simpmap,imname,logfile ; Display IDL header if help is required. if (keyword_set(help)) then begin get_idlhdr,'simpmap.pro' goto,finishup endif ; Make output file names. dotpos = STRPOS(imname,".") iflen = dotpos ifin1 = STRMID(imname,0,iflen) exlen = STRLEN(imname) - iflen ext = STRMID(imname,dotpos,exlen) hname = ifin1 + 'm' + ext ifin2 = STRMID(imname,1,iflen) sname = 'c' + ifin2 + ext ; Read input image im = READFITS(imname,imh,/NOSCALE) ; Find center, radius, p angle, b angle, etc. ew2 = FXPAR(imh,'CRVAL1') ns2 = FXPAR(imh,'CRVAL2') if ((ew2 gt 0) and (ns2 gt 0)) then dir = 'lr' if ((ew2 le 0) and (ns2 ge 0)) then dir = 'll' if ((ew2 lt 0) and (ns2 lt 0)) then dir = 'ul' if ((ew2 ge 0) and (ns2 le 0)) then dir = 'ur' datec = FXPAR(vih,'DATE-OBS') timec = FXPAR(vih,'TIME-OBS') READS(datec,day,mon,yrs,FORMAT='(i2,1x,i2,1x,i2)') READS(timec,hour,min,sec,FORMAT='(i2,1x,i2,1x,i2)') hr = hour + ((min + (sec / 60.0)) / 60.0) h = FXPAR(imh,'CENTERX') k = FXPAR(imh,'CENTERY') rad = FXPAR(imh,'RADIUS') h = h - 11 sunlong,mon,day,yrs,hr,p,b,lnaut,rsun,mjd dr = 1.74532925e-2 rd = 57.29577951 pr = p * dr br = b * dr lnr = lnaut * dr ; Scan image ; First make a subimage that only has the points we want to use. ims = im(10:503,*) ; Make arrays that contain the x and y distance from center Sun. xdist1 = (findgen(492) - h) ydist1 = (findgen(384) - k) xdsq = xdist1 ^ 2 ydsq = ydist1 ^ 2 ; Make an array that contains the radial distance from center Sun. ; and an array that contains the angle from north to east r = fltarr(492,384,/NOZERO) theta = fltarr(492,384,/NOZERO) for i = 0, 383 do begin r(0,i) = sqrt(xdsq + ydsq(i)) theta(0,i) = atan(xdist1, ydist1(i)) endfor negs = where(theta lt 0.0) theta(negs) = (2 * !PI) - theta(negs) ; free up memory no longer needed negs = 0 xdist1 = 0 ydist1 = 0 xdsq = 0 ydsq = 0 ; only want to work with pixels more than 15 arcsec in from limb. goodr = where(r le (rad-15*rad/s) latr = fltarr(492,384,/NOZERO) lonr = fltarr(492,384,/NOZERO) latr(goodr) = asin(sin(br) * cos(r(goodr)) $ + cos(br) * sin(r(goodr)) * cos(pr - theta(goodr)) lonr(goodr) = asin(sin(r(goodr)) $ * sin(pr - theta(goodr)) / cos(latr(goodr))) lonr = lonr + lnr lon = lonr * rd lat = latr * rd negs = where(lat lt 0.0) lat(negs) = -lat(negs) negs = where(lon lt 0.0) lon(negs) lon(negs) + 360.0 negs = 0 bigs = where(lon gt 360.0) lon(bigs) = lon(bigs) - 360.0 bigs = 0 ; Array operations cart = intarr(181,181) cart2 = intarr(181,181) cart3 = intarr(181,181) for i = 0, 491 do begin for j = 0, 384 do begin if (r(i,j) le (rad - 15*rad/s)) begin x = 90 + (lon - lnaut) y = 90 + lat cart(x,y) = cart(x,y) + im(i,j) cart2(x,y) = cart2(x,y) + 1 endif endfor endfor ;write logfile openw,lun,logfile, /GET_LUN for i = 0, 180 do begin for j = 0, 180 do begin if (cart2(i,j) ne 0) then begin cart3(i,j) = cart(i,j) / cart2(i,j) printf,lun,(i+lnaut-90),j-90,cart(i,j),cart2(i,j), $ cart3(i,j),FORMAT='(i4,i5,f15.4,f8.1,f13.4)' endif endfor endfor free_lun, lun ; Place arrays into images, place information in FITS headers, ; write FITS files. FXHMAKE,imouth,cart3 FXHMAKE,couth,cart2 FXWRITE,inamem,imouth,cart3 FXWRITE,cimname,couth,cart2 finishup: return end