;------------------------------------------------------------------ ;+ ; NAME: ; CENTELS ; 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. ; CATEGORY: ; POLAR ; CALLING SEQUENCE: ; centerls(vfile,jfile,imgscl,gcut) ; INPUTS: ; vfile = file name for input magnetogram (FITS format). ; jfile = file name for input direct image (FITS format). ; 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 in vfile and jfile 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 ;- ;------------------------------------------------------------------- pro centerls,vfinam,jfinam,imgscl,gcut,c_corr=gcorrk COMMON SUNRAD,sprad,quad ; Display IDL header if help is required. if (keyword_set(help)) then begin get_idlhdr,'centere.pro' goto,finishup endif ; Check for coordinate correction keyword if (keyword_set(gcorrk)) then gcorr=gcorrk else gcorr = 0.0 ; Read V file vi = READFITS(vfinam,vih,/NOSCALE) ji = READFITS(jfinam,jih,/NOSCALE) ; Obtain FITS parameters from V image ew = FXPAR(vih,'CRVAL1') ns = FXPAR(vih,'CRVAL2') 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) 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. ig = WHERE(xl1 GT 11 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 ; Write files CHECK_FITS,vi,vih,/UPDATE,/FITS CHECK_FITS,ji,jih,/UPDATE,/FITS WRITEFITS,vfinam,vi,vih WRITEFITS,jfinam,ji,jih 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