; ; WRITTEN: 15JUN2007 Xuepu Zhao for DRC ; MODIFIED: 24MAY2009 Xuepu Zhao ; pro ldrcmsc_updmsc20090706,crn,cl0,itl,updmsc0,updmscr,drcmsc,longi,xticscb,xticsct,xticsfb,xticsft,PEVT=pevt,DISP=disp,CRND=crnd,CLOG=clog if n_params( ) lt 1 then begin print,'ldrcmsc_updmsc20090706,crn,cl0,itl,updmsc0,updmscr,' print,' drcmsc,longi,...,/pevt' return end ; xsz=360 & ysz=180 & hxsz=xsz/4 sz=SIZE(updmsc0) & xsz=sz(1) & ysz=sz(2) & hxsz=xsz/4 ppd=xsz/360. ; Get sphd, slar2 and slar4 for calculating feature's differential rotation zgrid,xsz,ysz,phd,thd,lad,cth,sth sphd=strarr(xsz) for i=0,xsz-1 do begin CASE 1 OF phd(i) lt 10: sphd(i)='00'+strtrim(phd(i),2) phd(i) ge 10 and phd(i) lt 100: sphd(i)='0'+strtrim(phd(i),2) else: sphd(i)=strtrim(phd(i),2) ENDCASE endfor sphd=STRMID(sphd,0,7) slar2=fltarr(ysz) & slar4=slar2 for jj=0,ysz-1 do begin slarj=SIN(lad(jj)*!DTOR) slar2(jj)=slarj*slarj slar4(jj)=slarj^4. endfor ; Specify the reference time doy0 & hh0 (in days) for crn and cl0 scrn0=STRTRIM(crn,2) CASE 1 OF cl0 lt 10: scl0='00'+STRTRIM(cl0,2) cl0 ge 10 AND cl0 lt 100: scl0='0'+STRTRIM(cl0,2) else: scl0=STRTRIM(cl0,2) ENDCASE scrncl0=scrn0+':'+scl0 ; time_crncl0,scrncl0,symd_hms0 time_crncl0_ctimes,scrncl0,symd_hms0 year_sec_date,symd_hms0,yyyy0,mm0,dd0,hh0,mi0,ss0 dmy2doy,yyyy0,mm0,dd0,doy0 hh0=(hh0+mi0/60.+ss0/3600.)/24. fdoy0=doy0+hh0 ; Based on Pevtsov (2000) or Snograss (1983) formulae to get aa,bb,cc and ; the latitude of sunspot where synodic rate=27.2753 if keyword_set (pevt) then begin lat4carringtonrate,xsz,ysz,aa,bb,cc,lat2,ind,/pevt endif else begin lat4carringtonrate,xsz,ysz,aa,bb,cc,lat2,ind,/snog end ind0=ind(0) & ind1=ind(1) ; Get the Carrington longitude arraies sphdt0 & sphdtr for updmsc0 and updmscr dcl0itl=(cl0-itl)*ppd ;difference between cl0 and itl in pixels CASE 1 OF dcl0itl lt 0: begin sphd1l=sphd((dcl0itl+xsz):(xsz-1)) sphd1rl=STRTRIM(crn+1,2)+':'+sphd1l sphd2l=sphd(0:(xsz-1+dcl0itl)) sphd2rl=scrn0+':'+sphd2l sphdt0=[sphd1rl,sphd2rl] ;updmsc0 sphd1rl=scrn0+':'+sphd1l sphd2rl=STRTRIM(crn-1,2)+':'+sphd2l sphdtr=[sphd1rl,sphd2rl] ;updmscr end dcl0itl eq 0: begin sphdt0=scrn0+':'+sphd sphdtr=STRTRIM(crn-1,2)+':'+sphd end dcl0itl gt 0: begin sphd1l=sphd(dcl0itl:(xsz-1)) sphd1rl=scrn0+':'+sphd1l sphd2l=sphd(0:(dcl0itl-1)) sphd2rl=STRTRIM(crn-1,2)+':'+sphd2l sphdt0=[sphd1rl,sphd2rl] sphd1rl=STRTRIM(crn-1,2)+':'+sphd1l sphd2rl=STRTRIM(crn-2,2)+':'+sphd2l sphdtr=[sphd1rl,sphd2rl] end ENDCASE ; Correct the effect of differential rotation ; Get doyi,hhi and the time difference dt with erspect to doy0,hh0 longi=fltarr(2*xsz,ysz) & field=longi for ii=itl*ppd,xsz-1 do begin ; for ii=0,xsz-1 do begin scrncli=sphdt0(ii) ; time_crncl0,scrncli,symd_hmsi time_crncl0_ctimes,scrncli,symd_hmsi year_sec_date,symd_hmsi,yyyyi,mmi,ddi,hhi,mii,ssi dmy2doy,yyyyi,mmi,ddi,doyi hhi=(hhi+mii/60.+ssi/3600.)/24. fdoyi=doyi+hhi dt0=fdoyi-fdoy0 CASE 1 OF ABS(dt0) lt 55: dt=fdoyi-fdoy0 ; 0 > dt or dt > 0 in days dt0 ge 55: begin dr=year0 MOD 4 if dr eq 0 then dt=fdoyi-(dfoy0+366) else dt=fdoyi-(dfoy0+365) end else: begin dr=yeari MOD 4 if dr eq 0 then dt=(fdoyi+366)-dfoy0 else dt=(fdoyi+365)-fdoy0 end ENDCASE dlogCR1=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt ; dlogCR=(360/27.2753)*dt=13.1988*dt for jj=0,ysz-1 do begin dlogij10=(aa + bb*slar2(jj) + cc*slar4(jj))*dt dlogij1=dlogij10 - dlogCR1 ; longi(ii,jj)=phd(ii)-dlogij1 longi(ii,jj)=-dlogij10 ;phd(ii)=-dlogCR1 field(ii,jj)=updmsc0(ii,jj) endfor ;if ii gt xsz-5 then print,'dlogCR1,dlogij10,dlogij1:',dlogCR1,dlogij10,dlogij1 endfor is=ii ; Find the ending point for MSCR and get the contribution from MSCR dtt =fltarr(xsz) ; print,'sphdtr: ',sphdtr for ii=0,xsz-1 do begin scrncli=sphdtr(ii) ; time_crncl0,scrncli,symd_hmsi time_crncl0_ctimes,scrncli,symd_hmsi year_sec_date,symd_hmsi,yyyyi,mmi,ddi,hhi,mii,ssi dmy2doy,yyyyi,mmi,ddi,doyi hhi=(hhi+mii/60.+ssi/3600.)/24. fdoyi=doyi+hhi dt0=fdoyi-fdoy0 CASE 1 OF ABS(dt0) lt 55: dt=fdoyi-fdoy0 ; 0 > dt or dt > 0 in days dt0 ge 55: begin dr=year0 MOD 4 if dr eq 0 then dt=fdoyi-(dfoy0+366) else dt=fdoyi-(dfoy0+365) end else: begin dr=yeari MOD 4 if dr eq 0 then dt=(fdoyi+366)-dfoy0 else dt=(fdoyi+365)-fdoy0 end ENDCASE dtt(ii)=dt dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt dlogij20=(aa + bb*slar2(ysz-1) + cc*slar4(ysz-1))*dt dlogij2=dlogij20 - dlogCR ; longi_i179=360+phd(ii)-dlogij ; if longi_i179 gt 360 then goto,lbr if ii lt 5 then begin print,'dlogCR,dlogij20,dlogij2:',dlogCR,dlogij20,dlogij2 endif if ii gt xsz/2-5 then begin ; print,"scrncli,' ',symd_hmsi,dt0: ",scrncli,' ',symd_hmsi,dt0 print,'dlogCR,dlogij20,dlogij2:',dlogCR,dlogij20,dlogij2 endif if ABS(dlogij20) ge ABS(dlogCR1) then goto,lbr endfor lbr: iire=ii ;print,'dt:',dtt print,'is,iire:',is,iire ;return for ii=0,iire-1 do begin dt=dtt(ii) ; dlogCR=13.1988*dt dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt for jj=0,ysz-1 do begin dlogij0=(aa + bb*slar2(jj) + cc*slar4(jj))*dt dlogij=dlogij0 - dlogCR ; longi(ii+is,jj)=360+phd(ii)-dlogij longi(ii+is,jj)=-dlogij0 field(ii+is,jj)=updmscR(ii,jj) endfor endfor ; Get differential-rotation-corrected MSC, drcmsc drcmsc=fltarr(xsz+1,ysz) for j=0,ysz-1 do begin longi_midj0=longi(*,j) field_midj0=field(*,j) ind=WHERE(field_midj0 ne 0,cc) print,'cc:',cc longi_midj=longi_midj0(ind) field_midj=field_midj0(ind) for i=0,xsz do begin ; indi=where(longi_midj ge i-0.5 AND longi_midj lt i+0.5,ci) indi=where(longi_midj ge i-1.0 AND longi_midj lt i+1.0,ci) if ci gt 0 then begin br=TOTAL(field_midj(indi))/ci drcmsc(i,j)=br endif else begin drcmsc(i,j)=!values.f_nan endelse endfor indii=where(FINITE(drcmsc(*,j)) eq 1,cii) ; indin=where(FINITE(drcmsc(*,j)) ne 1,cin) ; if cii lt 361 then print,'j,cii:',j,cii if cii lt 361 then begin ti=indgen(361) xi=indii yi=drcmsc(indii,j) zi=SPLINE(xi,yi,ti) drcmsc(*,j)=zi endif endfor drcmsc=drcmsc(1:xsz,*) ; drcmscf='drcmscf_snog_'+'_'+STRTRIM(itl,2)+'_'+scrn0+'_'+STRTRIM(cl0,2)+'.fits' ; writefits,drcmscf,drcmsc,' ' end