;+ ; NAME : HR_CALIB.pro (procedure) ; PURPOSE : ; Calibrate the observed magnetic field data by the Huairou Solar Observation ; Station, correct for seeing (if set fseeing ne 1.0), show the vector magnetogram ; on the screen, or save it to a .ps file when set /saveps, or output the calibrated data ; to a FITS file given by FITSOUT keyword ; CATEGORY : ; idl/huairou ; CALLING SEQUENCE : ; hr_calib,lname,qname,uname,bx,by,bz,angle,head,nx=nx,ny=ny,sname=sname,$ ; saveps=saveps,dir=dir,ext=ext,nstep=nstep,magnify=magnify,xsize=xsize,$ ; fseeing=fseeing,ocenter=ocenter,qshift=qshift,ushift=ushift,/btonly,$ ; /old,/blonly,/btcont,/fitsout ; INPUTS : ; lname: file name of Stokes V profile data without extension ; qname: file name of Stokes Q profile data without extension ; uname: file name of Stokes U profile data without extension ; OPTIONAL INPUT PARAMETERS : ; none ; OUTPUTS : ; bx: tranverse magnetic field component in X direction (Gauss) ; by: tranverse magnetic field component in Y direction (Gauss) ; bz: longitudinal magnetic field (Gauss) ; angle: azimuthal angle of Bt ; head: head information. 12 by 3 or 4 (if set /sname) string array ; KEYWORD PARAMETERS : ; nx = nx: demanded pixel number in X direction. default: nx=170 ; ny = ny: demanded pixel number in Y direction. default: nx=120 ; old=old: if set, use the old method when calibrating Bt. ; i.e., Bt=60000*sqrt(Q^2+U^2). default: Bt=9370*sqrt(sqrt(Q^2+U^2)) ; dir=dir: directory where the observed data files are stored ; ext=ext: extension of the file names. default: '.dat' ; saveps : if set, save the aftermath in a PostScript file ; btonly : if set, show transversal field component only ; blonly : if set, show longitudinal field component only ; nstep : step increase to plot transversal field. default: nstep=10 ; sname : spot file name. if set, display spot image as background ; magnify: magnifying factor when showing the magnetogram. default: =3 ; btcont : if set, show the transversal field as contours ; fseeing: factor for seeing correct. default: =1.0 ; xsize : size (cm) in X direction when save PostScript file ; fitsout: if set, save the calibrated data in a FITS file. (head(*,0:2),bx,by,bz,angle) ; the file name is 'HByymmdd.hhmm.fit' where yymmdd.hhmm is the observing ; date and time in UT (hhmm) ; ocenter: two-element array containing the latitude & longitude of the center ; of the observed region (>0 for North and West). default: use the ; the value from file header provided by LNAME. ; qshift : two-element array containing required shift in X & Y direction ; for Q component. default = [0,0] ; ushift : two-element array containing required shift in X & Y direction ; for U component. default = [0,0] ; COMMON BLOCKS : none ; SIDE EFFECTS : ; RESTRICTIONS : ; PROCEDURE : ; ; MODIFICATION HISTORY : ; 2000/10/27 Hui LI, write ; 2001/05/21 Hui LI, add /ocenter, /qshift & /ushift ;------------------------------------------------------------------------ ;- pro hr_calib,lname,qname,uname,bx,by,bz,angle,head,nx=nx,ny=ny,old=old,sname=sname,$ saveps=saveps,dir=dir,ext=ext,btonly=btonly,blonly=blonly,nstep=nstep,$ magnify=magnify,btcont=btcont,fitsout=fitsout,fseeing=fseeing,xsize=xsize,$ ocenter=ocenter,qshift=qshift,ushift=ushift on_error,2 !except=2 if !version.os_family eq 'Windows' then olddev='win' $ else if !version.os_family eq 'unix' then olddev='x' else olddev=!d.name loadct,3 tvlct,0,0,0,!d.table_size-1 ; color of frame tvlct,0,0,0,!d.table_size-2 ; color of contour tvlct,255,255,255,0 ; back ground scale0=1 & flag=0 & scale=1.830 yp0=60 & x0=100 & y0=140 & imgmin=0.2 bfact=10. & plmin=0.002 c_labels=0 lname0=lname & qname0=qname & uname0=uname if keyword_set(sname) then sname0=sname ; ----------------- handle keywords ---------------------- if keyword_set(nx) then nx=nx else nx=170 if keyword_set(ny) then ny=ny else ny=120 if keyword_set(qshift) then qsh=qshift else qsh=[0,0] if keyword_set(ushift) then ush=ushift else ush=[0,0] if keyword_set(old) then begin cailt=3.0 method='OLD' endif else begin cailt=93.7/sqrt(2) method='NEW' endelse if keyword_set(ext) then ext=ext else ext='.dat' if keyword_set(nstep) then nstep=nstep else nstep=10 if keyword_set(fseeing) then fseeing=fseeing else fseeing=1.0 if keyword_set(magnify) then magnify=magnify else $ if !d.x_size gt 400 then magnify=3 else magnify=2.5 if keyword_set(dir) then begin if !version.os_family eq 'Windows' then begin if strmid(dir,strlen(dir)-1,1) eq '\' then dir=dir else dir=dir+'\' endif else begin if strmid(dir,strlen(dir)-1,1) eq '/' then dir=dir else dir=dir+'/' endelse endif else begin dir='' endelse ; ------------ Determine and read files ---------------------- lname=dir+lname+ext qname=dir+qname+ext uname=dir+uname+ext ff=findfile(lname) !error=0 if ff(0) eq '' then begin !error=1 if keyword_set(dir) then $ print,'File '+lname+' not found ! Procedure stopped.' else $ print,'File '+lname+' in current directory not found ! Procedure stopped.' return endif ff=findfile(qname) !error=0 if ff(0) eq '' then begin !error=1 if keyword_set(dir) then $ print,'File '+qname+' not found ! Procedure stopped.' else $ print,'File '+qname+' in current directory not found ! Procedure stopped.' return endif ff=findfile(uname) !error=0 if ff(0) eq '' then begin !error=1 if keyword_set(dir) then $ print,'File '+uname+' not found ! Procedure stopped.' else $ print,'File '+uname+' in current directory not found ! Procedure stopped.' return endif hr_read,lname,ll0,hd1 hr_read,qname,qq0,hd2 hr_read,uname,uu0,hd3 qq0=shift(qq0,qsh(0),qsh(1)) uu0=shift(uu0,ush(0),ush(1)) ;print,'I detect the center of the observed area is located at:' ;print,' latitude : '+hd1(4)+' (>0 for north; <0 for south)' ;print,' longitude: '+hd1(5)+' (>0 for west; <0 for east)' ;print,'Do you want to correct the location ? (Y/N) (N)' ;anw=get_kbrd(1) ;if strupcase(anw) eq 'Y' then begin ; print,'Please input the latitude and longitude :' ; read,latit,longit ; delX=longit-float(hd1(5)) ; delY=latit-float(hd1(4)) ; hd1(4)=string(latit,format='(f5.1)') ; hd1(5)=string(longit,format='(f5.1)') ; hd1(8)=string(float(hd1(8))+dely,format='(f5.1)') ; hd1(9)=string(float(hd1(9))+delx,format='(f5.1)') ;endif if keyword_set(saveps) or keyword_set(fitsout) then begin tmp=strupcase(strmid(hd1(1),0,3)) if tmp eq 'JAN' then mon='01' if tmp eq 'FEB' then mon='02' if tmp eq 'MAR' then mon='03' if tmp eq 'APR' then mon='04' if tmp eq 'MAY' then mon='05' if tmp eq 'JUN' then mon='06' if tmp eq 'JUL' then mon='07' if tmp eq 'AUG' then mon='08' if tmp eq 'SEP' then mon='09' if tmp eq 'OCT' then mon='10' if tmp eq 'NOV' then mon='11' if tmp eq 'DEC' then mon='12' date=strmid(hd1(1),9,2)+mon+strmid(hd1(1),4,2) if keyword_set(old) then tmp='HOB' else tmp='HNB' ofname=dir+tmp+date+'.'+strmid(hd1(2),0,2)+strmid(hd1(2),3,2) endif if keyword_set(saveps) then begin scale0=25 & flag=1 ps_file=ofname+'.ps' ss=size(saveps) if ss(ss(0)+1) eq 7 then begin sl=strlen(saveps)-2 if (sl gt 1) and (strupcase(strmid(saveps,sl,2)) eq 'PS') then $ ps_file=dir+saveps else ps_file=dir+saveps+'.ps' endif endif if keyword_set(ocenter) then begin latit=ocenter(0) & longit=ocenter(1) delX=longit-float(hd1(5)) delY=latit-float(hd1(4)) hd1(4)=string(latit,format='(f5.1)') hd1(5)=string(longit,format='(f5.1)') hd1(8)=string(float(hd1(8))+dely,format='(f5.1)') hd1(9)=string(float(hd1(9))+delx,format='(f5.1)') endif head=[[hd1],[hd2],[hd3]] gdat=fltarr(nx,ny) qq=gdat & uu=qq for i=0,nx-1 do begin for j=0,ny-1 do begin gdat(i,j)=mean(ll0(3*i:3*i+2,4*j:4*j+3)) qq(i,j)=mean(qq0(3*i:3*i+2,4*j:4*j+3)) uu(i,j)=mean(uu0(3*i:3*i+2,4*j:4*j+3)) endfor endfor if keyword_set(sname) then begin sname=dir+sname+ext ff=findfile(sname) !error=0 if ff(0) eq '' then begin !error=1 if keyword_set(dir) then $ print,'File '+sname+' not found ! Procedure stopped.' else $ print,'File '+sname+' in current directory not found ! Procedure stopped.' return endif hr_read,sname,img0,hd4 head=[[hd1],[hd2],[hd3],[hd4]] img=congrid(float(img0),nx,ny,/inter,/min) endif ; -------------------- Calculate longitudinal field strength ----------------- x=[20.,40.,80.,160.,280.,424.,556.,677.,788.,890.,984.,$ 1071.,1152.,1320.,1500.,1700.,2000.] y=[20.,40.,80.,160.,320.,640.,960.,1280.,1600.,1920.,2240.,$ 2560.,2880.,3335.,3665.,3900.,4000.] x=congrid(x,100,/inter,/min) ;added by Hui LI on 2000/10/25 y=smooth(congrid(y,100,/inter,/min),3) ;added by Hui LI on 2000/10/25 caill=0.5 gdat=gdat*caill bz=fltarr(nx,ny) bz=interpol(y,x,gdat) ii=where(gdat lt 0.0,count) if count ne 0 then bz(ii)=-interpol(y,x,-gdat(ii)) bz=reform(bz,nx,ny)*fseeing bzmax=max(bz) & bzmin=min(bz) ; ------------- Calculate transversal field strength and direction ----------- angle=fltarr(nx,ny) rdat=sqrt(qq^2+uu^2) if keyword_set(old) then bt=cailt*rdat $ else bt=cailt*sqrt(rdat) btmax=max(bt*fseeing) ii=where( qq eq 0. and uu ge 0. , count ) if count ne 0 then angle(ii)=!pi/2./2 ii=where( qq eq 0. and uu lt 0. , count ) if count ne 0 then angle(ii)=3*!pi/2./2 ii=where( qq ne 0. , count ) if count ne 0 then angle(ii)=!pi/2.-0.5*atan(uu(ii),qq(ii)) bx=bt*cos(angle)*fseeing by=bt*sin(angle)*fseeing ; >> determine the Bt direction << if nx eq 1 or ny eq 1 then return print,'getting potential field ...' cff1n,congrid(bz,128,128,/inter,/min),bxp,byp bxp=congrid(bxp,nx,ny,/inter,/min) byp=congrid(byp,nx,ny,/inter,/min) print,'setting Bt direction...' flip = where( bx*bxp + by*byp lt 0.0, count ) if count gt 0 then begin bx( flip ) = -bx( flip ) by( flip ) = -by( flip ) angle( flip ) = angle(flip)+!pi endif angle=angle*180./!pi print,'Bt_max = ',btmax,' (Gauss)' print,'Bz_max = ',bzmax,' (Gauss) Bz_min = ',bzmin,' (Gauss)' ;-------------- Open window ------------------- ;if not keyword_set(fitsout) then begin if keyword_set(neutral) then begin if keyword_set(nocaption) then $ if keyword_set(noplot) then offset=-80 else offset=20 $ else $ if keyword_set(noplot) then offset=-30 else offset=120 endif else begin if keyword_set(nocaption) then offset=-100 else offset=0 endelse y0=offset+yp0 wwx=nx*magnify+2*x0 ;expected window size wwy=ny*magnify+y0+170 if !d.window ne -1 then wno=!d.window else wno=0 if not keyword_set(printout) and not keyword_set(saveps) then begin if (!d.window eq -1) or (wwx gt !d.x_size+1) or (wwy gt !d.y_size+1) then $ window,wno,xsize=wwx,ysize=wwy pos=[scale0*x0,scale0*y0,scale0*(x0+nx*magnify-1),scale0*(y0+ny*magnify-1)] endif else begin if keyword_set(xsize) then xsize=float(xsize) else xsize=18.0 ysize=xsize*ny/nx pos=[0,0,fix(xsize*1000),fix(ysize*1000)] xoff=fix((21-xsize)/2.+0.5) & yoff=fix((29.7-ysize)/2+0.5) set_plot,'ps' device,file=ps_file,xs=xsize,ys=ysize,xoff=xoff,yoff=yoff endelse erase ; ------------------- Display Spot -------------------------- if keyword_set(sname) then begin imgmax=max(img) if not keyword_set(nospot) then begin img=sqrt(sqrt(img)) img=(img-0.9*min(img))/max(img-0.9*min(img))*240 if keyword_set(saveps) then begin img1=congrid(img,2*nx,2*ny,/interp) tv,img1,pos(0),pos(1),/dev,xsize=pos(2)-pos(0),ysize=pos(3)-pos(1) endif else begin img1=congrid(img,pos(2)-pos(0),pos(3)-pos(1),/interp) tv,img1,pos(0),pos(1),/dev endelse endif endif ;------ display contour of longitudinal or transversal fields ------- if not keyword_set(btonly) then begin if keyword_set(btcont) then begin r=smooth(congrid(bt,nx,ny,/interp),2) levp=[200,500,800,1200,1600,2000,2500,3000] levstr='level=' endif else begin r=float(smooth(congrid(bz,nx,ny,/interp),2)) levp=[40,80,160,320,640,960,1280,1600,1920,2240,2560,2880] levstr='level=!9+!3' endelse nn=n_elements(levp) for i=0,nn-2 do $ levstr=levstr+strcompress(string(levp(i),format='(i)'),/remove_all)+',' levstr=levstr+string(levp(nn-1),format='(I4)')+' Gauss' if not keyword_set(btcont) then begin levn=levp for i=0,nn-1 do levn(i)=-levp(nn-i-1) lev=[levn,levp] lnum=n_elements(lev) lstyle=intarr(lnum) lthick=intarr(lnum) lstyle(where(lev lt 0))=1 ; dot lthick(where(lev lt 0))=2 if not keyword_set(sname) then begin lcolor=replicate(!d.n_colors/1.9,lnum) lcolor(where(lev lt 0))=!d.n_colors/1.7 ;(where(lev lt 0))=300 endif else lcolor=replicate(!d.n_colors-10,lnum) endif else begin lev=levp lstyle=intarr(n_elements(lev)) lthick=intarr(n_elements(lev)) endelse lstyle(where(lev gt 0))=0 ; solid line lthick(where(lev gt 0))=1 if not keyword_set(nolong) then begin if flag eq 1 then $ contour,r,pos=pos,/dev,xstyle=1+4,ystyle=1+4,xrange=[0.,nx],$ yrange=[0.,ny],level=lev,c_linestyle=lstyle,c_thick=lthick,$ /noerase,spline=spline,c_labels=c_labels,back=255,color=1, $ xtitle=xtitle,ytitle=ytitle $ else $ contour,r,pos=pos,/dev,xstyle=1+4,ystyle=1+4,xrange=[0.,nx],$ yrange=[0.,ny],level=lev,c_linestyle=lstyle,c_thick=lthick,$ xtitle=xtitle,ytitle=ytitle, $ /noerase,spline=spline,c_labels=c_labels,back=255,c_color=lcolor ; use follow method and surpress labels endif endif ;---------------- display transversal fields -------------- if keyword_set(old) then btmin=10 else btmin=30 if keyword_set(saveps) then clr=10 else $ if keyword_set(sname) then clr=!d.n_colors/7 else $ if !d.n_colors lt 255 then clr=10 else clr=!d.n_colors-100 if not keyword_set(blonly) then begin if flag eq 1 then plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase,color=clr, $ xrange=[0.,nx*magnify],yrange=[0.,ny*magnify],xstyle=1+4,ystyle=1+4 else $ plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase,color=150, $ xrange=[0.,nx*magnify],yrange=[0.,ny*magnify],xstyle=1+4,ystyle=1+4 nstep=max([nstep/magnify,1]) ;initially set nstep=10 for j=nstep,ny-nstep,nstep do begin for i=nstep,nx-nstep,nstep do begin ; if img(i,j) gt imgmax*imgmin then begin Bt1=(bt(i,j)>1)/bfact plen=sqrt((Bt1-Btmin)>1) ;Bfact=10,Btmin=100 if Bt1 gt Btmin then begin th=atan(by(i,j),bx(i,j)) x=plen*cos(th)/magnify y=plen*sin(th)/magnify if !order eq 0 then begin x1=[i-x,i+x] y1=[j-y,j+y] endif else begin x1=[i-x,i+x] y1=[ny-j-1+y,ny-j-1-y] endelse oplot,x1*magnify,y1*magnify,color=clr dth=20./180.*!pi len=sqrt((x1(1)-x1(0))^2+(y1(1)-y1(0))^2) larrow=min([len/4.,5.]) ddx=larrow*(x1(0)-x1(1))/len ddy=larrow*(y1(0)-y1(1))/len xx1=x1(1)+ddx*cos(dth)+ddy*sin(dth) yy1=y1(1)-ddx*sin(dth)+ddy*cos(dth) xx2=x1(1)+ddx*cos(-dth)+ddy*sin(-dth) yy2=y1(1)-ddx*sin(-dth)+ddy*cos(-dth) oplot,[xx1,x1(1),xx2,xx1]*magnify, $ [yy1,y1(1),yy2,yy1]*magnify,color=clr endif ;endif endfor endfor endif ;--------------- Plot frame ------------------- title='SOLAR VECTOR MAGNETOGRAM, BEIJING' xtitle='arc sec' & ytitle='arc sec' xmax=314 & ymax=218 ;Field of view: 5.23 by 3.63 arc min minxy=min([xmax,ymax]) lg=fix(alog10(minxy)) step=10.^lg tickv=fltarr(10) for i=0,9 do tickv(i)=step*i xrange=[0.,xmax] & yrange=[0.,ymax] datime='Date & Time: '+hd1(1)+' '+hd1(2) if keyword_set(blonly) then begin files='Files Name: '+strmid(lname,0,8) cont='Contents : Longitudinal Field' endif else begin if keyword_set(btonly) then begin files='Files Names: '+strmid(qname,0,8)+' '+strmid(uname,0,8) cont='Contents : Transversal Field' endif else begin files='Files Names: '+strmid(lname,0,8)+' '+$ strmid(qname,0,8)+' '+strmid(uname,0,8) cont='Contents : Vector Magnetogram' endelse endelse tmp1=strcompress(hd1(4),/re) & tmp2=strcompress(hd1(5),/re) if float(hd1(4) gt 0) then sn='N'+strcompress(hd1(4),/remove) $ else sn='S'+strmid(tmp1,1,strlen(tmp1)-1) if float(hd1(5) gt 0) then ew='W'+strcompress(hd1(5),/remove) $ else ew='E'+strmid(tmp2,1,strlen(tmp2)-1) sunc='Sun Disk Coord. : '+sn+' '+ew orth='Orthogonal Coor. : '+hd1(6)+' '+hd1(7) carr='Carrington Coor. : '+hd1(8)+' '+hd1(9) wave='Wave Length : '+strcompress(hd1(0),/remove) see='Seeing : '+strcompress(hd1(3),/remove) fram='Frame Number : '+strcompress(hd1(10),/remove) if keyword_set(saveps) then begin plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase, $ xrange=xrange,yrange=yrange,xstyle=1,ystyle=1, $ xticks=fix(nx*scale/step),xtickv=tickv,xminor=5,xticklen=0.03, $ yticks=fix(ny*scale/step),ytickv=tickv,yminor=5,yticklen=0.03, $ xtitle=xtitle,ytitle=ytitle,charsize=charsize xyouts,pos(0)*0.5,pos(3)+4500,/dev,title,chars=2.5,color=3 xyouts,pos(0)+1000,pos(3)+3500,/dev,datime,chars=1.,color=3 xyouts,pos(0)+1000,pos(3)+3000,/dev,files,chars=1,color=3 xyouts,pos(0)+1000,pos(3)+2500,/dev,cont,chars=1,color=3 xyouts,pos(0)+1000,pos(3)+2000,/dev,sunc,chars=1,color=3 xyouts,pos(0)+1000,pos(3)+1500,/dev,orth,chars=1,color=3 xyouts,pos(0)+1000,pos(3)+1000,/dev,carr,chars=1,color=3 xyouts,pos(2)*0.6,pos(3)+3500,/dev,wave,chars=1,color=3 xyouts,pos(2)*0.6,pos(3)+3000,/dev,see,chars=1,color=3 xyouts,pos(2)*0.6,pos(3)+2500,/dev,fram,chars=1,color=3 if not keyword_set(btonly) then $ xyouts,pos(0)+1000,pos(3)+500,/dev,levstr,chars=1,color=3 xyouts,1000,-yoff*850,/dev,'Created by IDL at '+string(systime()),chars=1,color=3 device,/close print,'Result has been saved in PostScript file: '+ps_file set_plot,olddev endif else begin plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase, $ xrange=xrange,yrange=yrange,xstyle=1,ystyle=1, $ xticks=fix(nx*scale/step),xtickv=tickv,xminor=5,xticklen=0.03, $ yticks=fix(ny*scale/step),ytickv=tickv,yminor=5,yticklen=0.03, $ xtitle=xtitle,ytitle=ytitle,color=clr,charsize=charsize xyouts,pos(0)*0.5,pos(3)+130,/dev,title,chars=2.,color=clr xyouts,pos(0),pos(3)+105,/dev,datime,chars=0.8,color=clr xyouts,pos(0),pos(3)+90,/dev,files,chars=0.8,color=clr xyouts,pos(0),pos(3)+75,/dev,cont,chars=0.8,color=clr xyouts,pos(0),pos(3)+60,/dev,sunc,chars=0.8,color=clr xyouts,pos(0),pos(3)+45,/dev,orth,chars=0.8,color=clr xyouts,pos(0),pos(3)+30,/dev,carr,chars=0.8,color=clr xyouts,pos(2)*0.7,pos(3)+105,/dev,wave,chars=0.8,color=clr xyouts,pos(2)*0.7,pos(3)+90,/dev,see,chars=0.8,color=clr xyouts,pos(2)*0.7,pos(3)+75,/dev,fram,chars=0.8,color=clr if not keyword_set(btonly) then $ xyouts,pos(0),pos(3)+15,/dev,levstr,chars=0.8,color=clr endelse ;endif ; ----------- Output FITS file ----------------------- if keyword_set(fitsout) then begin if keyword_set(sname) then array=fltarr(nx,ny,4) $ else array=fltarr(nx,ny,3) array(*,*,0)=bx array(*,*,1)=by array(*,*,2)=bz if keyword_set(sname) then array(*,*,3)=congrid(float(img0),nx,ny,/inter,/min) mkhdr,hdr,array sxaddpar,hdr,'LNAME',lname0,'File name of L-component',/pdu sxaddpar,hdr,'QNAME',qname0,'File name of Q-component',/pdu sxaddpar,hdr,'UNAME',uname0,'File name of U-component',/pdu if keyword_set(sname) then $ sxaddpar,hdr,'SNAME',sname0,'File name of the spot image',/pdu sxaddpar,hdr,'WAVELENT',hd1(0),'Observation wavelenght',/pdu sxaddpar,hdr,'DATE-OBS',hd1(1),'Date of observation',/pdu sxaddpar,hdr,'TIME-OBS',hd1(2),'Time of observation in UT',/pdu sxaddpar,hdr,'SEEING',hd1(3),'Seeing on the Day of observation',/pdu sxaddpar,hdr,'INSTRUME','Solar Vector Magnetograph','Observation Instrument',/pdu sxaddpar,hdr,'OBSERVER','Huairou, Beijing Astro. Obs.','Observer',/pdu sxaddpar,hdr,'DISK-COX',float(hd1(5)),'Sun disk coordinate X',/pdu sxaddpar,hdr,'DISK-COY',float(hd1(4)),'Sun disk coordinate Y',/pdu sxaddpar,hdr,'ORTH-COX',float(hd1(6)),'Orthogonal coordinate X',/pdu sxaddpar,hdr,'ORTH-COY',float(hd1(7)),'Orthogonal coordinate Y',/pdu sxaddpar,hdr,'CARR-COX',float(hd1(9)),'Carrington coordinate X',/pdu sxaddpar,hdr,'CARR-COY',float(hd1(8)),'Carrington coordinate Y',/pdu sxaddpar,hdr,'FRAMEs',fix(hd1(10)),'Integration frames',/pdu sxaddpar,hdr,'HR-NUM',fix(hd1(11)),'Huairou sunspot number',/pdu sxaddpar,hdr,'SEE-FACT',fseeing,'Seeing correction factor. 1.0 for uncorrected',/pdu sxaddpar,hdr,'METHOD',method,'Calibrating method: OLD or NEW',/pdu sxaddpar,hdr,'PROJECT',0,' 0/1: uncorrected/corrected for projection effect',/pdu writefits,ofname+'.fit',array,hdr endif ;---------- end