; CODE FROM JESPER SCHOU ; SLIGHTLY MODIFIED BY S. COUVIDAT ;------------------------------------------------------------------------------------------------ PRO cubicfit,X,A,F,pder F = A[0]*X+A[1]*X^2.d0+A[2]*X^3.d0 pder= [ [X],[X^2.d0],[X^3.d0] ] END PRO linearfit,X,A,F,pder F = A*X pder= [X] END PRO seb_sun_lin,FSN1 ; FSN1 is the first FSN of the non-linearity sequence ; the last FSN is FSN1+55 ; SEQUENCES TAKEN FROM SPACE ground=0 ; will only work on sequences taken from space ;SPAWN,'show_info ds="hmi.lev1[][36793811-36793866]" -qp > temp' ;SPAWN,'show_info ds="hmi.lev1[][36793811-36793866]" -q key=HCAMID,EXPTIME > temp2' ;SPAWN,'show_info ds="hmi.lev1[][36782291-36782346]" -qp > temp' ;SPAWN,'show_info ds="hmi.lev1[][36782291-36782346]" -q key=HCAMID,EXPTIME > temp2' ;SPAWN,'show_info ds="hmi.lev1[][36785171-36785226]" -qp > temp' ;SPAWN,'show_info ds="hmi.lev1[][36785171-36785226]" -q key=HCAMID,EXPTIME > temp2' ;SPAWN,'show_info ds="hmi.lev1[][4649267-4649322]" -qp > temp' ;SPAWN,'show_info ds="hmi.lev1[][4649267-4649322]" -q key=HCAMID,EXPTIME > temp2' ;SPAWN,'show_info ds="hmi.lev1[][4649411-4649466]" -qp > temp' ;SPAWN,'show_info ds="hmi.lev1[][4649411-4649466]" -q key=HCAMID,EXPTIME > temp2' ;SPAWN,'show_info ds="hmi.lev1_nrt[][62818424-62818479]" -qp > temp' ; HFTSACID=3034 ;SPAWN,'show_info ds="hmi.lev1_nrt[][62818424-62818479]" -q key=HCAMID,EXPTIME > temp2' ;SPAWN,'show_info ds="hmi.lev1_nrt[][62821304-62821359]" -qp > temp' ; HFTSACID=3034 ;SPAWN,'show_info ds="hmi.lev1_nrt[][62821304-62821359]" -q key=HCAMID,EXPTIME > temp2' ; SEQUENCES TAKEN FROM THE GROUND ;SPAWN,'show_info ds="hmi_ground.lev0[1420880-1420945]" -qp seg=file > temp' ;SPAWN,'show_info ds="hmi_ground.lev0[1420880-1420945]" -q key=HMI_SEQ_ID_EXP_PATH,HMI_FSW_IMG_CMDED_EXPOSURE > temp2' ;ground=1 SPAWN,'show_info ds="hmi.lev1[]['+STRTRIM(STRING(FSN1),1)+'-'+STRTRIM(STRING(FSN1+55),1)+']" -qp > temp' ; HFTSACID=3034 SPAWN,'show_info ds="hmi.lev1[]['+STRTRIM(STRING(FSN1),1)+'-'+STRTRIM(STRING(FSN1+55),1)+']" -q key=HCAMID,EXPTIME > temp2' SPAWN,'show_info ds="hmi.lev1[]['+STRTRIM(STRING(FSN1),1)+']" key=T_OBS -q',TOBS PRINT,'BEGINNING TIME: ',TOBS IF(ground EQ 1) THEN nfiles=66 ELSE nfiles=56 nbin=64 filenames=STRARR(nfiles) data=FLTARR(2,nfiles) OPENR,1,'temp' READF,1,filenames CLOSE,1 OPENR,1,'temp2' READF,1,data CLOSE,1 if (ground EQ 1) THEN data[1,*]=data[1,*]/1000. ; convert exposure times from milliseconds to seconds image=FLTARR(nbin,nbin,nfiles) IF( ground EQ 1) THEN BEGIN darkf=FLTARR(4096,4096) darks=darkf FOR i=0,2 DO darkf+=FITSIO_READ_IMAGE(filenames[i]) FOR i=0,2 DO darks+=FITSIO_READ_IMAGE(filenames[i+33]) darkf=darkf/3. darks=darks/3. darkf=REBIN(darkf,nbin,nbin) darks=REBIN(darks,nbin,nbin) ENDIF PRINT,'READING FILES' FOR i=0,nfiles-1 DO BEGIN PRINT,i if (ground EQ 1) THEN d=FITSIO_READ_IMAGE(filenames[i]) ELSE d=FITSIO_READ_IMAGE(filenames[i]+'/image_lev1.fits') if (ground EQ 1) THEN image[*,*,i]=REBIN(d,nbin,nbin) ELSE image[*,*,i]=REBIN(d,nbin,nbin)*data[1,i] if (ground EQ 1 AND data[0,i] EQ 1) THEN image[*,*,i]-=darkf ;remove dark frame on gound data if (ground EQ 1 AND data[0,i] EQ 0) THEN image[*,*,i]-=darks ENDFOR if (ground EQ 1) THEN front=WHERE(data[0,*] EQ 1) ELSE front=WHERE(data[0,*] EQ 3) if (ground EQ 1) THEN side=WHERE(data[0,*] EQ 0) ELSE side=WHERE(data[0,*] EQ 2) expf=data[1,front]*1000. exps=data[1,side]*1000. f00=image[31,31,front] f01=image[31,32,front] f10=image[32,31,front] f11=image[32,32,front] fxx=(f00+f01+f10+f11)/4.0 ; Find brightest part intmax=0 for i=0,nbin-1 do for j=0,nbin-1 do intmax=intmax>image[i,j,front]>image[i,j,side] af=WHERE(intmax LE 12000,naf) ;pfxx=poly_fit(expf[af],fxx[af],1) ; linear intensity for front camera W=FLTARR(naf)+1.0 A=50. res=curvefit(expf[af],fxx[af],W,A,FUNCTION_NAME='linearfit') pfxx=[0.,A] PRINT,pfxx s00=image[31,31,side] s01=image[31,32,side] s10=image[32,31,side] s11=image[32,32,side] sxx=(s00+s01+s10+s11)/4.0 as=WHERE(intmax LE 12000,nas) W=FLTARR(nas)+1.0 A=50. res=curvefit(exps[as],sxx[as],W,A,FUNCTION_NAME='linearfit') psxx=[0.,A] PRINT,psxx W=FLTARR(naf)+1.0 A=[-20.,3.d-3,-1.d-7] res=curvefit(fxx[af],(fxx[af]-poly(expf[af],pfxx))*1000.d0,W,A,FUNCTION_NAME='cubicfit',ITER=iter,TOL=1.d-12,STATUS=status,ITMAX=200) resf=[0.,A[0]/1000.d0,A[1]/1000.d0,A[2]/1000.d0] PRINT,resf,iter,status W=FLTARR(nas)+1.0 A=[-20.,3.d-3,-1.d-7] res=curvefit(sxx[as],(sxx[as]-poly(exps[as],psxx))*1000.d0,W,A,FUNCTION_NAME='cubicfit',ITER=iter,TOl=1.d-12,STATUS=status,ITMAX=200) ress=[0.,A[0]/1000.d0,A[1]/1000.d0,A[2]/1000.d0] PRINT,ress,iter,status PRINT,'PLOTTING RESULTS' SET_PLOT,'ps' DEVICE,FILE='non_linearity_'+STRTRIM(STRING(FSN1),1)+'_updated.ps',xoffset=0,yoffset=0,xsize=21,ysize=22,/color LOADCT,3 !P.MULTI=[0,1,2] IF (GROUND EQ 0) THEN PLOT,expf,fxx/1000.,psym=2,thick=3,tit='!17 TOBS='+TOBS,xtit='exposure time (ms)',ytit='actual average int. (10!U3!N DN)',charsize=1.5,yst=1,xst=1,xrange=[0,300],yrange=[0,17] else PLOT,expf,fxx/1000.,psym=2,thick=3,tit='!17 TOBS='+TOBS,xtit='exposure time (ms)',ytit='actual average int. (10!U3!N DN)',charsize=1.5,yst=1,xst=1,xrange=[0,17000],yrange=[0,17] OPLOT,FINDGEN(3000000)/100.,POLY(FINDGEN(3000000)/100.,pfxx)/1000.,col=180,thick=2 OPLOT,exps,sxx/1000.,psym=4,thick=3 OPLOT,FINDGEN(3000000)/100.,POLY(FINDGEN(3000000)/100.,psxx)/1000.,col=180,thick=2,linestyle=2 legend,['front camera','side camera'],psym=[2,4] yo='front:'+STRJOIN(STRING(REFORM(pfxx[1])),/SINGLE) xyouts,20,12,yo yo='side:'+STRJOIN(STRING(REFORM(psxx[1])),/SINGLE) xyouts,20,11,yo ;IF (GROUND EQ 1) THEN plot,fxx/1000.,fxx-poly(expf,pfxx)-offsetf,xst=1,thick=3,psym=2,tit='!17 Intensity averaged over entire CCD',xtit='actual intensity (10!U3!N DN)',ytit='actual intensity - linear intensity (DN)',charsize=1.5,yst=1,xrange=[0,15],yrange=[-50,150] ELSE plot,fxx/1000.,fxx-poly(expf,pfxx)-offsetf,xst=1,thick=3,psym=2,tit='!17 Intensity averaged over entire CCD',xtit='actual intensity (10!U3!N DN)',ytit='actual intensity - linear intensity (DN)',charsize=1.5,yst=1,xrange=[0,15],yrange=[-100,100] IF (GROUND EQ 1) THEN plot,fxx/1000.,fxx-poly(expf,pfxx),xst=1,thick=3,psym=2,tit='!17'+TOBS,xtit='actual intensity (10!U3!N DN)',ytit='actual intensity - linear intensity (DN)',charsize=1.5,yst=1,xrange=[0,15],yrange=[-50,150] ELSE plot,fxx/1000.,fxx-poly(expf,pfxx),xst=1,thick=3,psym=2,tit='!17',xtit='actual intensity (10!U3!N DN)',ytit='actual int. - linear int. (DN)',charsize=1.5,yst=1,xrange=[0,15],yrange=[-50,150] legend,['front camera','side camera'],psym=[2,4] ;oplot,sxx/1000.,sxx-poly(exps,psxx)-offsets,psym=4,thick=3 oplot,sxx/1000.,sxx-poly(exps,psxx),psym=4,thick=3 oplot,FINDGEN(15000)/1000.,POLY(FINDGEN(15000),ress),col=180,thick=2,linestyle=2 oplot,FINDGEN(15000)/1000.,POLY(FINDGEN(15000),resf),col=180,thick=2 yo='front:'+STRJOIN(STRING(REFORM(resf[1:3])),/SINGLE) IF(GROUND EQ 0) THEN xyouts,0.2,100,yo ELSE xyouts,0.2,100,yo yo='side:'+STRJOIN(STRING(REFORM(ress[1:3])),/SINGLE) IF(GROUND EQ 0) THEN xyouts,0.2,90,yo ELSE xyouts,0.2,90,yo OPLOT,[12,12],[-1000,1000],thick=2 DEVICE,/CLOSE PRINT,'non-linearity coeffs for front camera=',resf PRINT,'non-linearity coeffs for side camera=',ress SET_PLOT,'X' !P.MULTI=0 END