;
; WRITTEN:  15JUN2007   Xuepu Zhao for DRC 
; MODIFIED: 24MAY2009   Xuepu Zhao
;
pro Ldrcmsc_updmsc,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_updmsc,crn,cl0,itl,updmsc0,updmscR,drcmsc,longi,'
       print,'              /pevt,/disp'
       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
    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
; Set the reference time doy0 & hh0 (in day) for the reference time
    scrn0=STRTRIM(crn,2)
    scl0=STRTRIM(cl0,2)
    scrncl0=scrn0+':'+scl0    
    time_crncl0,scrncl0,symd_hms0
    yyyy0=FIX(strmid(symd_hms0,0,4))
    mm0=FIX(strmid(symd_hms0,5,2))
    dd0=FIX(strmid(symd_hms0,8,2))
    dmy2doy,yyyy0,mm0,dd0,doy0
    hh0=(FIX(strmid(symd_hms0,11,2))+FIX(strmid(symd_hms0,15,2))/60.)/24.
; Based on Pevtsov (2000) or Snofrass (1983) formulae to get the latitude
; of sunspot that is used by Carrington to obtain CR rate of 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 longitude difference from cl0
    sphd=strtrim(phd,2)
    dcl0itl=(cl0-itl)*ppd
    CASE 1 OF
     dcl0itl lt 0: begin
      phd1=phd((dcl0itl+xsz):(xsz-1))
      sphd1=STRTRIM(crn+1,2)+':'+STRTRIM(phd1,2)
      phd2=phd(0:(xsz-1+dcl0itl))
      sphd2=scrn0+':'+STRTRIM(phd2,2)
      sphdt0=[sphd1,sphd2]
      sphd1=scrn0+':'+STRTRIM(phd1,2)
      sphd2=STRTRIM(crn-1,2)+':'+STRTRIM(phd2,2)
      sphdtr=[sphd1,sphd2]
                   end
     dcl0itl eq 0: begin
      sphdt0=scrn0+':'+sphd
      sphdtr=STRTRIM(crn-1,2)+':'+sphd
                   end
     dcl0itl gt 0: begin
      phd1=phd(dcl0itl:(xsz-1))
      sphd1=scrn0+':'+STRTRIM(phd1,2)
      phd2=phd(0:(dcl0itl-1))
      sphd2=STRTRIM(crn-1,2)+':'+STRTRIM(phd2)
      sphdt0=[sphd1,sphd2]
      sphd1=STRTRIM(crn-1,2)+':'+STRTRIM(phd1,2)
      sphd2=STRTRIM(crn-2,2)+':'+STRTRIM(phd2,2)
      sphdtr=[sphd1,sphd2]
                   end
    ENDCASE
; Correct the effect of differential rotation
; Get doyi,hhi and the time difference dt with erspect to the reference time
    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
        yyyyi=FIX(strmid(symd_hmsi,0,4))
        dmy2doy,yyyyi,FIX(strmid(symd_hmsi,5,2)),$
                FIX(strmid(symd_hmsi,8,2)),doyi
        hhi=FIX(strmid(symd_hmsi,11,2))/24.
        dt=(doyi-doy0)+(hhi-hh0)     ; 0 > dt > 0 in days
        for jj=0,ysz-1 do begin
            dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
        ;   dlogCR=13.1988*dt
            dlogij=(aa + bb*slar2(jj) + cc*slar4(jj))*dt - dlogCR
            longi(ii,jj)=phd(ii)-dlogij
            field(ii,jj)=updmsc0(ii,jj)
        endfor
    endfor
    is=ii+1
; Find the ending point for MSCR and get the contribution from MSCR
    dtt =fltarr(xsz)
    for ii=0,xsz-1 do begin
        scrncli=sphdtr(ii)
        time_crncl0,scrncli,symd_hmsi
        yyyyi=FIX(strmid(symd_hmsi,0,4))
        dmy2doy,yyyyi,FIX(strmid(symd_hmsi,5,2)),$
                FIX(strmid(symd_hmsi,8,2)),doyi
        hhi=FIX(strmid(symd_hmsi,11,2))/24.      ; in day
        dt=(doyi-doy0)+(hhi-hh0)                 ; dt < 0 in days for MSCR
        dtt(ii)=dt
        dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
        dlogij=(aa + bb*slar2(ysz-1) + cc*slar4(ysz-1))*dt - dlogCR 
        longi_i179=360+phd(ii)-dlogij
        if longi_i179 gt 360 then goto,lbr
    endfor
  lbr: iire=ii 
    for ii=0,iire-1 do begin
        dt=dtt(ii)
        for jj=0,ysz-1 do begin
            dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
        ;   dlogCR=13.1988*dt
            dlogij=(aa + bb*slar2(jj) + cc*slar4(jj))*dt - dlogCR
            longi(ii+is,jj)=360+phd(ii)-dlogij
            field(ii+is,jj)=updmscR(ii,jj)
        endfor
    endfor
  ; xszr=xsz+iire-1
  ; xe=360+phd(ii) & ie=is
; Get differential-rotation-corrected MSC, drcmsc
    drcmsc=fltarr(xsz+1,ysz)
    for j=0,ysz-1 do begin
        longi_midj=longi(*,j)
        field_midj=field(*,j)
        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='Ldrcmscf'+'_'+scrn0+'_'+STRTRIM(cl0,2)+'.fits'
    writefits,drcmscf,drcmsc,' '
; Display
    if keyword_set(disp) then begin
       mx=100
       tek_color
       outfile='drcmsc_'+scrncl0+'.new.ps'
       psout,041,outfile
       !p.multi=[0,1,4]
    ;  ynon = REPLICATE(' ', 2)
    ;  ytic = 
       xnon = REPLICATE(' ', 2)
       xtic = ['1937:290',xnon,'1936:20',xnon,'1936:110',xnon,'1936:200',xnon,'1936:290',xnon,'1935:20',xnon,'1935:110']
       xn=n_elements(xtic)
       xtn=xn-1
       dx=360/xtn               ;Instead of 540 see pscimg06.pro
       cxtv = FINDGEN(xn)*dx
       txnon= replicate(' ',xn)
; Panel 1: Updmscs 
       org3msc=longi
       org3msc(0:89,*)=updmscL(270:359,*)
       org3msc(90:449,*)=updmsc0
       org3msc(450:539,*)=updmscR(0:89,*)
       borg3msc=bytscl(org3msc,max=mx,min=-mx)
       pscimg06,borg3msc,'Wide MDI Synoptic Carrington Map Centered at CT 1636:200',px,py,sy='sla'
       AXIS,XAXIS=0,XTICKS=xtn,XTICKV=cxtv,XTICKN=xtic,xtit='Carrington (blue) Longitude',xcharsize=1
       AXIS,XAXIS=1,XTICKS=xtn,XTICKV=cxtv,XTICKN=txnon
       coe=360./540
       for i=0,539,10 do oplot,[i,i]*coe,[-1,1],color=4
     ; oplot,[cl0,cl0],[-1,1],thick=5,color=2
       oplot,[270,270]*coe,[-1,1],thick=5,color=2
       i0=90 & i360=449 
       oplot,[i0,i0]*coe,[-1,1],thick=3
       oplot,[i360,i360]*coe,[-1,1],thick=3
;
       lad=SIN(lad*!DTOR)            ; 12May2006
       bbn0=fltarr(xszr,180)
       bbn=bbn0
       bbn(xszl:xsz0-1,*)=drcmsc
       bdrcmsc=bytscl(bbn,max=mx,min=-mx)
;Panel 2
       dx=540/xtn               ;Instead of 540 see pscimg06.pro
       cxtv = FINDGEN(xn)*dx-90
       txnon= replicate(' ',xn)
       plot,longi(0,*),lad,/nodata,psym=3,yrange=[-1,1],xstyle=5,$
       ;    ytitle='Heliolatitude',xrange=[0,540],ystyle=1,$
            ytitle='Sine Heliolatitude',xrange=[-90,450],ystyle=1,$
       ;    xtitle='Hybrid (Carrington and Hleographic) longitude', $
            title="Stretching of Carrington Longitudes (blue lines) to Heliographic Locations (Colored Curves)"
       AXIS,XAXIS=0,XTICKS=xtn,XTICKV=cxtv,XTICKN=xtic,xtit='Carrington (blue) Longitude',xcharsize=1
       AXIS,XAXIS=1,XTICKS=xtn,XTICKV=cxtv,XTICKN=txnon
       for ii=0,xszl-1,10 do oplot,longi(ii,*),lad,psym=3,color=2
       for ii=xszl-1,xsz0-1,10 do oplot,longi(ii,*),lad,psym=3
       for ii=xsz0+9,xszr-1,10 do oplot,longi(ii,*),lad,psym=3,color=2
       for i=-90,449,10 do oplot,[i,i],[-1,1],color=4
       oplot,[-90,449],[sinl0,sinl0]
       oplot,[-90,449],[sinl1,sinl1] 
; Panel 3
       xtic = ['98.05.23:290',xnon,'98.05.23:20',xnon,'98.05.23:110',xnon,'98.05.23:200',xnon,'98.05.23:290',xnon,'98.05.23:20',xnon,'98.05.23:110']
       dx=540/xtn               ;Instead of 540 see pscimg06.pro
       cxtv = FINDGEN(xn)*dx-90
       txnon= replicate(' ',xn)
       clong=fltarr(1.5*xsz,ysz)
       for ii=-90,449,10 do begin
           dhli=(ii-180)            
           for jj=0,ysz-1 do begin
               dt=dhli/(aa + bb*slar2(jj) + cc*slar4(jj))
               dlogCR=13.1988*dt
               clong(ii+90,jj)=180+dlogCR
           endfor
       endfor
       plot,clong(0,*),lad,/nodata,psym=3,yrange=[-1,1],xstyle=5,$
       ;    ytitle='Heliolatitude',xrange=[0,540],ystyle=1,$
            ytitle='Sine Heliolatitude',xrange=[-90,450],ystyle=1,$
       ;    xtitle='Hybrid (Carrington and Hleographic) longitude', $
            title="Pulling from Carrington Locations (Blue Curves) to Heliographic Longitude (Black Vertical Lines)"
       AXIS,XAXIS=0,XTICKS=xtn,XTICKV=cxtv,XTICKN=xtic,xtit='Heliographic (black) Longitude',xcharsize=1
       AXIS,XAXIS=1,XTICKS=xtn,XTICKV=cxtv,XTICKN=txnon
       for i=-90,449,10 do oplot,[i,i],[-1,1]
     ; for ii=0,xszl-1,10 do oplot,longi(ii,*),lad,psym=3,color=4
       for ii=-90,449,10 do oplot,clong(ii+90,*),lad,psym=3,color=4
     ; for ii=xsz0+9,xszr-1,10 do oplot,longi(ii,*),lad,psym=3,color=4
       oplot,[-90,449],[sinl0,sinl0]
       oplot,[-90,449],[sinl1,sinl1]
; Panel 4
       xtic = ['98.05.23:290',xnon,'98.05.23:20',xnon,'98.05.23:110',xnon,'98.05.23:200',xnon,'98.05.23:290',xnon,'98.05.23:20',xnon,'98.05.23:110']
       dx=360/xtn               ;Instead of 540 see pscimg06.pro
       cxtv = FINDGEN(xn)*dx
       txnon= replicate(' ',xn)
       bbn0=fltarr(xszr,180)
       bbn=bbn0
       bbn(xszl:xsz0-1,*)=drcmsc
       bdrcmsc=bytscl(bbn,max=mx,min=-mx)
       pscimg06,bdrcmsc,'MDI Synchronic Map for 1998.05.23_16:03:30',$
             px,py,sy='sla'
       AXIS,XAXIS=0,XTICKS=xtn,XTICKV=cxtv,XTICKN=xtic,xtit='Heliographic (black) Longitude',xcharsize=1
       AXIS,XAXIS=1,XTICKS=xtn,XTICKV=cxtv,XTICKN=txnon
       oplot,[cl0,cl0],[-90,90],thick=3,color=2
       for i=90,449,10 do oplot,[i,i]*coe,[-1,1]
     ; sc_annot,1,xtic=xtic,longi='heliographic',/norm
     ; sc_annot,1,longi='heliographic',/norm
       !p.multi=0
       psdone
       print,outfile
     endif
end