PRO DISKCOR3,sphe,sphcor,xdiam,ydiam,shft,xc,yc ; ; Version 3 - without the spectroheliogram rotation ; and with the optional output of shft,xc,yc ; ;Gives a circular shape to the spectroheliogram, i.e. the correction of the ; horizontal diameter of the disk and the compensation for the entrance slit curvature. ; INPUT: sphe - input spectroheliogram ; OUTPUT: sphcor - corrected spectroheliogram (without annotation) ; xdiam,ydiam - x,y diameters of the solar disk ; shft - an array of shifts in x to rectify the disk ; xc,yc - coordinates of the disk centre ; ; Copyright M. Klvana and M. Sobotka, ASU AVCR Ondrejov, Czech Rep., 2009 ; si=size(sphe) sph=sphe ; threshold to define the disk window,0 beg=fix(total(sph)/n_elements(sph)/10) ;beg = min distance of secondary maximum in bins, ;calculated as the mean intensity histy,sph,ii,hist,bin=10,/abs ma=max(hist,pos1) ;1st max. near zero ma=max(hist(beg:*),pos2) pos2=pos2+beg mi=min(hist(pos1:pos2),mipos) ;minimum between the two maxima mipos=mipos+pos1 thresh=ii(mipos) plots,[ii(pos1),ii(pos1)],[0,ma],col=0 plots,[ii(pos2),ii(pos2)],[0,ma],col=0 plots,[ii(mipos),ii(mipos)],[0,ma],col=0 print,'Disk edge intensity:',thresh ; find the centroid of disk mask=sph gt thresh w=where(mask,co) ;1-D coordinates x=w mod si(1) ;2-D coords y=w/si(1) xc=total(x)/co ;x-coordinate of centroid yc=total(y)/co ;y-coordinate of centroid ; x and y diameters of the disk beltx=total(sph(*,yc-50:yc+50),2)/101. beltx=beltx < ii(pos2) beltx=smooth(beltx,3) belty=total(sph(xc-50:xc+50,*),1)/101. belty=belty < ii(pos2) belty=smooth(belty,3) beltdx=beltx-shift(beltx,1) beltdx=beltdx(1:*) ma=max(beltdx,pa) mi=min(beltdx,pi) xdiam=abs(pi-pa) ;x-diameter of solar disk beltdy=belty-shift(belty,1) beltdy=beltdy(1:*) ma=max(beltdy,pa) mi=min(beltdy,pi) ydiam=abs(pi-pa) ;y-diameter of solar disk aspect=float(ydiam)/float(xdiam) ;current aspect ratio print,'X and Y diameters:',xdiam,ydiam ; making the Sun roundish - correction of the x-diameter sph=congrid(sph,round(si(1)*aspect),si(2),cubic=(-0.5)) si=size(sph) ; new size ; correction for entrance slit curvature mask=sph gt thresh shfts=fltarr(si(2)) for j=0,si(2)-1 do begin w=where(mask(*,j),co) if co gt 0 then shfts(j)=(w(co-1)+w(0))/2. endfor ws=where(shfts) ; parabolic fit through shifts yy=findgen(si(2)) parfit=poly_fit(yy(ws),shfts(ws),2) shft=parfit(1)*yy+parfit(2)*yy^2 ; spectroheliogram rectification sph2=sph for j=0,si(2)-1 do begin xx=findgen(si(1))+shft(j) rect=interpolate(sph(*,j),xx) sph2(*,j)=rect endfor ; find the centroid of rectified disk mask=sph2 gt thresh w=where(mask,co) ;1-D coordinates x=w mod si(1) ;2-D coords y=w/si(1) xc=long(total(x)/co) ;x-coordinate of centroid yc=long(total(y)/co) ;y-coordinate of centroid ; shift the centroid to the centre of image sphcor=shift(sph2,si(1)/2-xc,si(2)/2-yc) ; visualization window,0,xs=si(1),ys=si(2) sph3=sphcor ang=dindgen(361)*2*!pi/360 x=cos(ang)*(ydiam+6)/2+si(1)/2 y=sin(ang)*(ydiam+6)/2+si(2)/2 sph3(x,y)=max(sphcor) tvscl,sph3 wait,2 END