;------------------------------------------------------------------ ;+ ; NAME: ; CENLSVAR ; PURPOSE: ; Find center Sun position in BBSO polar magnetograms. ; This procedure estimates the Sun center position in pixel ; coordinates. This procedure uses the center X, Y position ; (provided by the guider) along with image scale information ; provided by the user. ; ; This version allows for the fact that the number of bad ; columns on the left side of a J image may be variable! ; CATEGORY: ; POLAR ; CALLING SEQUENCE: ; cenlsvar(vimg,jimg,vih,jih,imgscl,gcut) ; INPUTS: ; vi = IDL image for input magnetogram ; ji = IDL image for input direct image ; vih = FITS header for vimg ; jih = FITS header for jimg ; imgscl = image scale (pixels / arcsec). ; gcut = gradient cutoff value ; gcorr = correction to picture center coordinates. The ; real coordinates are the file coordinates divided ; by gcorr. If gcorr is zero or not present, no correction ; is applied. ; KEYWORD PARAMETERS: ; OUTPUTS: ; The FITS headers vih and jih have the following ; FITS keywords added: ; CENTERX, CENTERY = center Sun coordinates in pixels. ; RADIUS = radius of Sun in pixels. ; SDCX = standard deviation of CENTERX in pixels ; SDCY = standard deviation of CENTERY in pixels ; ; COMMON BLOCKS: ; SUNRAD contains: ; R = radius of Sun in pixels. ; Q = quadrant of circle to fit ; NOTES: ; The time, date, and position information are taken from the ; V image. The FITS keywords are added in BOTH images. ; This assumes that the time and position are the same for both ; images. ; MODIFICATION HISTORY: ; J. Varsik, 18 Mar. 1996 ; J. Varsik, 12 Jan. 1999 -- numerous changes. ;- ;------------------------------------------------------------------- pro cenlsvar,vi,ji,vih,jih,imgscl,gcut,c_corr=gcorrk COMMON SUNRAD,sprad,quad ; Display IDL header if help is required. if (keyword_set(help)) then begin get_idlhdr,'cenlsvar.pro' goto,finishup endif ; Check for coordinate correction keyword if (keyword_set(gcorrk)) then gcorr=gcorrk else gcorr = 0.0 ; Obtain FITS parameters from V image ew = FXPAR(vih,'CRVAL1') ns = FXPAR(vih,'CRVAL2') datec = FXPAR(vih,'DATE-OBS') timec = FXPAR(vih,'TIME-OBS') IF (STRMID(datec, 2, 1) EQ '/') THEN BEGIN READS,datec,day,mon,yrs,FORMAT='(i2,1x,i2,1x,i2)' yrs = yrs + 1900 ENDIF ELSE BEGIN READS,datec,yrs,mon,day,FORMAT='(i4,1x,i2,1x,i2)' ENDELSE READS,timec,hour,min,sec,FORMAT='(i2,1x,i2,1x,i2)' hr = hour + ((min + (sec / 60.0)) / 60.0) ; Use solar ephemeris from Johns Hopkins library sun,yrs,mon,day,hr,pa=p,lat0=b,sd=rsun ; pandb,mon,day,yrs,hr,p,b,rsun,mjd sprad = rsun * imgscl sdrad = 0.0 if (gcorr EQ 0.0) then begin h = imgscl * ew + 256 k = (-1.0 * imgscl * ns) + 192 endif else begin h = imgscl * (ew / gcorr) + 256 k = (-1.0 * imgscl * (ns / gcorr)) + 192 endelse print,h,k ; now we have our first estimate of the Sun center position. ; First we get the quadrant for the Sun center. This is passed ; to the function nslimb in a common block. quad = 0 IF (h LT 0.0) AND (abs(k) LT abs(h)) THEN quad = 2 IF (h GT 0.0) AND (abs(k) LT abs(h)) THEN quad = 4 IF (k LT 0.0) AND (abs(h) LE abs(k)) THEN quad = 1 IF (k GE 0.0) AND (abs(h) LE abs(k)) THEN quad = 3 IF quad EQ 0 THEN begin print,'Error assigning quadrant' return ENDIF ; Next we isolate the limb points. We use the Sobel function to ; give the gradient of the direct image. jgrad = SOBEL(ji) ; Next we use the region-coloring function to identify all the ; regions of the direct image with large gradients. jgcut = jgrad GT gcut jglabel = LABEL_REGION(jgcut) jgmag = HISTOGRAM(jglabel,REVERSE_INDICES=r) ; Now we select the largest region. This should be the set of limb ; points. imsiz = SIZE(ji) xax = imsiz(1) jgmagmax = MAX(jgmag(1:*),maxreg) maxreg = maxreg + 1 p = r(r(maxreg):r(maxreg+1)-1) xl1 = INTARR(N_ELEMENTS(p)) yl1 = INTARR(N_ELEMENTS(p)) FOR i = 0, N_ELEMENTS(p)-1 DO BEGIN kpix = p(i) xl1(i) = kpix MOD xax yl1(i) = kpix / xax ENDFOR ; Be sure to avoid bad points. ; If north or south limb, adjust left bad column badcol = 11 IF (quad eq 1) THEN BEGIN ; north limb badrow = jglabel(0:200,30) bcs = where (badrow eq maxreg) badcol = max(bcs) ENDIF IF (quad eq 3) THEN BEGIN ; south limb badrow = jglabel(0:200,364) bcs = where (badrow eq maxreg) badcol = max(bcs) ENDIF ig = WHERE(xl1 GT badcol AND xl1 LT 504) xl = xl1(ig) yl = yl1(ig) ; Now we use the non-linear least squares curve-fitting function ; to locate the new Sun center points. a = FLTARR(2) a(0) = h a(1) = k w = MAKE_ARRAY(N_ELEMENTS(ig),/FLOAT,VALUE=1.0) IF (quad EQ 1) OR (quad EQ 3) THEN BEGIN yfit = CURVEFIT(xl,yl,w,a,siga,FUNCTION_NAME='NSLIMB') ENDIF ELSE BEGIN xfit = CURVEFIT(yl,xl,w,a,siga,FUNCTION_NAME='EWLIMB') ENDELSE h = a(0) k = a(1) sigh = siga(0) sigk = siga(1) PRINT,h,sigh,k,sigk ; New we put the center Sun coordinates in the FITS headers. commx = 'center of sun, x-coordinate' commy = 'center of sun, y-coordinate' commr = 'radius of sun, pixels' commr2 = 'standard deviation in CENTERX' commr3 = 'standard deviation in CENTERY' FXADDPAR,vih,'CENTERX',h,commx FXADDPAR,vih,'CENTERY',k,commy FXADDPAR,vih,'RADIUS',sprad,commr FXADDPAR,vih,'SDCX',sigh,commr2 FXADDPAR,vih,'SDCY',sigk,commr3 FXADDPAR,jih,'CENTERX',h,commx FXADDPAR,jih,'CENTERY',k,commy FXADDPAR,jih,'RADIUS',sprad,commr FXADDPAR,jih,'SDCX',sigh,commr2 FXADDPAR,jih,'SDCY',sigk,commr3 finishup: RETURN END PRO NSLIMB,x,a,f,pder COMMON SUNRAD,r,q ; function for use with CURVEFIT. ; r is set to solar radius in pixels. CASE q OF 1: BEGIN f = a(1) + SQRT(r^2 - (x - a(0))^2) IF N_PARAMS() GE 4 THEN BEGIN pder = FLTARR(N_ELEMENTS(x),2) pder(*,0) = (x - a(0)) / SQRT(r^2 - (x - a(0))^2) pder(*,1) = 1.0 ENDIF END 3: BEGIN f = a(1) - SQRT(r^2 - (x - a(0))^2) IF N_PARAMS() GE 4 THEN BEGIN pder = FLTARR(N_ELEMENTS(x),2) pder(*,0) = -1.0*(x - a(0)) / SQRT(r^2 - (x - a(0))^2) pder(*,1) = 1.0 ENDIF END ENDCASE END PRO EWLIMB,y,a,f,pder COMMON SUNRAD,r,q ; function for use with CURVEFIT. ; r is set to solar radius in pixels. CASE q OF 2: BEGIN f = a(0) + SQRT(r^2 - (y - a(1))^2) IF N_PARAMS() GE 4 THEN BEGIN pder = FLTARR(N_ELEMENTS(y),2) pder(*,0) = 1.0 pder(*,1) = (y - a(1)) / SQRT(r^2 - (y - a(1))^2) ENDIF END 4: BEGIN f = a(0) - SQRT(r^2 - (y - a(1))^2) IF N_PARAMS() GE 4 THEN BEGIN pder = FLTARR(N_ELEMENTS(x),2) pder(*,0) = 1.0 pder(*,1) = -1.0*(y - a(1)) / SQRT(r^2 - (y - a(1))^2) ENDIF END ENDCASE END