CAP update

Waypoint in pipeline update

c2jtimeline.pro

@@ -14,7 +14,7 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
   ;
   ; IDL> c2jtimeline, /pickpath
   ;
-  ; c2jtimeline'2 detailed results are plotted into .ps files, one devoted to the EXLIN results, one to the EXCUB results.
+  ; c2jtimeline's detailed results are plotted into .ps and jpg files, one devoted to the EXLIN results, one to the EXCUB results.
   ; On screen, by default only the whole dataset/timeline is shown. To also display on screen the detailed plots,
   ; each showing the cnt2Jy trends inside a single time interval, use:
   ;
@@ -25,7 +25,7 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
   ; orders of magnitude of 10E+06, 10E+07, producing properly-formatted output tables.
   ;
   ; Authors: Marcello Giroletti, Simona Righini
-  ; Last edited: Nov 3, 2017
+  ; Last edited: Jan 20, 2022
   ;
 
   sep=path_sep()
@@ -117,8 +117,17 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
 
     if f eq 0 then origin='EXLIN' else origin='EXCUB'
 
-    real0=where(cnt2Jy_0_el ne -99)
-    real1=where(cnt2Jy_1_el ne -99)
+    real0=where(cnt2Jy_0_el gt 0)
+    real1=where(cnt2Jy_1_el gt 0)
+    nreal0=n_elements(real0)
+    nreal1=n_elements(real1)
+    c2j0stats=moment(cnt2Jy_0_el[real0])
+    c2j1stats=moment(cnt2Jy_1_el[real1])
+    c2j0av=c2j0stats[0]
+    c2j1av=c2j1stats[0]
+    c2j0stdev=sqrt(c2j0stats[1])
+    c2j1stdev=sqrt(c2j1stats[1])
+
     wholeplot=errorplot(el_d[real0], cnt2Jy_0_el[real0], el_d[real0]*0, err_cnt2Jy_0_el[real0], name='CH0')
     wholeplot.title='All cnt2Jy measurements - '+origin
     wholeplot.xtitle='Elevation (deg)'
@@ -135,8 +144,29 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
     wholeplot2.symbol='triangle'
     wholeplot2.sym_filled=1
     wholeplot2.sym_size=1.0
+    wholeplot3=plot([min(el_d[real0]),max(el_d[real0])], [c2j0av,c2j0av],/overplot, name='Aver0')
+    wholeplot3.linestyle='-'
+    wholeplot3.color='midnight blue'
+    wholeplot4=plot([min(el_d[real1]),max(el_d[real1])], [c2j1av,c2j1av],/overplot, name='Aver1')
+    wholeplot4.linestyle='-'
+    wholeplot4.color='deep sky blue'
+    wholeplot5=plot([min(el_d[real0]),max(el_d[real0])], [c2j0av+c2j0stdev,c2j0av+c2j0stdev],/overplot, name='StDev0')
+    wholeplot5.linestyle='dash_dot'
+    wholeplot5.color='midnight blue'
+    wholeplot6=plot([min(el_d[real0]),max(el_d[real0])], [c2j0av-c2j0stdev,c2j0av-c2j0stdev],/overplot)
+    wholeplot6.linestyle='dash_dot'
+    wholeplot6.color='midnight blue'
+    wholeplot7=plot([min(el_d[real1]),max(el_d[real1])], [c2j1av+c2j1stdev,c2j1av+c2j1stdev],/overplot, name='StDev1')
+    wholeplot7.linestyle='dash_dot'
+    wholeplot7.color='deep sky blue'
+    wholeplot8=plot([min(el_d[real1]),max(el_d[real1])], [c2j1av-c2j1stdev,c2j1av-c2j1stdev],/overplot)
+    wholeplot8.linestyle='dash_dot'
+    wholeplot8.color='deep sky blue'
+
 
     leg = legend(TARGET=[wholeplot, wholeplot2], POSITION=[0.9,0.9], /AUTO_TEXT_COLOR)
+
+    wholeplot.save, workpath+'timeline_'+origin+'.jpg'
 
 
     ; overall plot of full dataset
@@ -259,7 +289,7 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
           ;          print, ' '
           ;          stop
 
-          realc2j0=where(c2J0plot ne -99)
+          realc2j0=where(c2J0plot gt 0)
           ystep=(max(c2J0plot[realc2j0]+errc2J0plot[realc2j0])-min(c2J0plot[realc2j0]-errc2J0plot[realc2j0]))/10.0
           labstep=ystep
           if ystep eq 0 then begin
@@ -293,8 +323,8 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
               displot.title=origin+' cnt2Jy_0 vs Time, INT'+strcompress(string(j),/remove_all)
               displot.xtitle='Elapsed time (hh.h)'
               displot.ytitle='Jy/cnt'
-              displot.xrange=[min(th0plot)-xstep,max(th0plot)+xstep]
-              displot.yrange=[min(c2J0plot[realc2j0])-ystep,max(th0plot[realc2j0])+ystep]
+           ;   displot.xrange=[min(th0plot)-xstep,max(th0plot)+xstep]
+           ;   displot.yrange=[min(c2J0plot[realc2j0])-ystep,max(th0plot[realc2j0])+ystep]
             endif
           endif
 
@@ -392,7 +422,7 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
           y_1min=Ch_1[0]+(Ch_1[1])*th_1val+Ch_sig1[0]
           y_1max=Ch_1[0]+(Ch_1[1])*th_1val-Ch_sig1[0]
 
-          realc2j1=where(c2J1plot ne -99)
+          realc2j1=where(c2J1plot gt 0)
           ystep=(max(c2J1plot[realc2j1]+errc2J1plot[realc2j1])-min(c2J1plot[realc2j1]-errc2J1plot[realc2j1]))/10.0
           labstep=ystep
           if ystep eq 0 then begin
@@ -426,8 +456,8 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
               displot6.title=origin+' cnt2Jy_1 vs Time, INT'+strcompress(string(j),/remove_all)
               displot6.xtitle='Elapsed time (hh.h)'
               displot6.ytitle='Jy/cnt'
-              displot6.xrange=[min(th1plot)-xstep,max(th1plot)+xstep]
-              displot6.yrange=[min(c2J1plot[realc2j1])-ystep,max(th1plot[realc2j1])+ystep]
+            ;  displot6.xrange=[min(th1plot)-xstep,max(th1plot)+xstep]
+            ;  displot6.yrange=[min(c2J1plot[realc2j1])-ystep,max(th1plot[realc2j1])+ystep]
             endif
           endif
 
@@ -482,6 +512,8 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
         case mode of
           'aver': begin
             printf, Unit, FORMAT = '(i4," ",d8.1,2(d11.4),8(g13.7))',j+1,freq[0],ti,tf,0,x_0mean,0,sigma_0m,0,x_1mean,0,sigma_1m
+            print, x_0mean, sigma_0m, sigma_0d
+            print, x_1mean, sigma_1m, sigma_1d
           end
           'linf': begin
             printf, Unit, FORMAT = '(i4," ",d8.1,2(d11.4),8(g13.7))',j+1,freq[0],ti,tf,C_0[1],C_0[0],C_sig0[1],C_sig0[0],C_1[1],C_1[0],C_sig1[1],C_sig1[0]
@@ -500,10 +532,10 @@ pro c2jtimeline, pickpath=pickpath, range=range, linf=linf, detplot=detplot
   device, /close
   close, /ALL
 
-  print, ' ******************************'
-  print, ' ****  C2JTIMELINE is DONE ****'
-  print, ' *** Next step is RUNTARGET ***'
-  print, ' ******************************'
+  print, ' ****************************** '
+  print, ' ****  C2JTIMELINE is DONE **** '
+  print, ' *** Next step is RUNTARGET *** '
+  print, ' ****************************** '
 
   return
 end

calibfit.pro

@@ -122,6 +122,7 @@ pro calibfit, scanflag,stacflag,polyflag,section,tipo,allpath,namefile,Out3,flux
       rescut=[residual[0:resrange],residual[-1*resrange,-1]]   ; avoiding the central part of the subscan, where artifacts can be present due to sidelobes
       resstat=moment(rescut)
       res_rms=sqrt(resstat[1])
+     ; err_cnt=sqrt(err_cnt^2+res_rms^2)   ; XXX temporaneo 
       err_cnt=sqrt(err_cnt^2+res_rms^2+(0.03*peak_cnt)^2)  ; updated error for the amplitude measurement, including a default 3% uncertainty on calibrator flux-amplitude
 
       getcnt2Jy, flux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, Out3, Level
@@ -157,6 +158,7 @@ pro calibfit, scanflag,stacflag,polyflag,section,tipo,allpath,namefile,Out3,flux
       rescut_b=[residual[0:resrange],residual[-1*resrange,-1]]   ; avoiding the central part of the subscan, where artifacts can be present due to sidelobes
       resstat=moment(rescut_b)
       res_rms=sqrt(resstat[1])
+     ; err_cnt=sqrt(err_cnt^2+res_rms^2)  ; XXX eliminazione temporanea 3% sull'errore del calibratore
       err_cnt=sqrt(err_cnt^2+res_rms^2+(0.03*peak_cnt)^2)  ; updated error for the amplitude measurement
 
       getcnt2Jy, flux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, Out3, Level

runcalib.pro

@@ -158,7 +158,7 @@ pro cal_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
   gaintime=data[0].time              ; associated MJD (first sample of first subscan)
   firstscanname=strsplit(sublist[0],sep,/extract)
   firstscandate=strsplit(firstscanname[-1],'-',/extract)
-  
+
   datesplit=strsplit(firstscandate[0],path_sep(),/extract)
   yyyymmdd=datesplit[-1]
 
@@ -1226,12 +1226,16 @@ pro cal_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         peak_cnt_0 = cal_wmean(p0, e0, /nan)
         peak_cnt_1 = cal_wmean(p1, e1, /nan)
-        w0=1/(e0/p0)
-        w1=1/(e1/p1)
-        err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
-        err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
-        ;  err_cnt_0 = stddev(p0, /nan)
-        ;  err_cnt_1 = stddev(p1, /nan)
+        w0=1.0/(e0/p0)
+        w1=1.0/(e1/p1)
+        ;  err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
+        ;  err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
+        err_cnt_0 = 1.0/sqrt((1.0/e0[0])^2+(1.0/e0[1])^2)   ; aggiornati 29/11/2021
+        err_cnt_1 = 1.0/sqrt((1.0/e1[0])^2+(1.0/e1[1])^2)
+      ; PROVA: 
+      ;  err_cntsd_0 = stddev(p0, /nan)
+      ;  err_cntsd_1 = stddev(p1, /nan)
+      ;  print, 'Ehilà: ', err_cnt_0, err_cntsd_0
 
         ; riscalo i ct2Jy e gli errori della polarizzazione 0 per l'offset incrociato
         c0[0]=c0[0]*exp(-(offset[1]*1.66/beamd)^2.)
@@ -1246,11 +1250,12 @@ pro cal_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         cnt2Jy_0 = cal_wmean(c0, d0, /nan)
         cnt2Jy_1 = cal_wmean(c1, d1, /nan)
-        w0=1/(d0/c0)
-        w1=1/(d1/c1)
-        err_cnt2Jy_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1])
-        err_cnt2Jy_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1])
-
+        w0=1.0/(d0/c0)
+        w1=1.0/(d1/c1)
+        ; err_cnt2Jy_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1])
+        ; err_cnt2Jy_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1])
+        err_cnt2Jy_0 = 1.0/sqrt((1.0/d0[0])^2+(1.0/d0[1])^2)   ; aggiornati 29/11/2021
+        err_cnt2Jy_1 = 1.0/sqrt((1.0/d1[0])^2+(1.0/d1[1])^2)
 
 
         ; last check: when calibrators are not known, the resulting dummy cnt2Jy values are to be reset to -99.00
@@ -1326,10 +1331,12 @@ pro cal_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         peak_cnt_0 = cal_wmean(p0c, e0c, /nan)
         peak_cnt_1 = cal_wmean(p1c, e1c, /nan)
-        w0=1/(e0c/p0c)
-        w1=1/(e1c/p1c)
-        err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
-        err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
+        w0=1.0/(e0c/p0c)
+        w1=1.0/(e1c/p1c)
+        ;  err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/sqrt(w0[0]+w0[1])
+        ;  err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/sqrt(w1[0]+w1[1])
+        err_cnt_0 = 1.0/sqrt((1.0/e0c[0])^2+(1.0/e0c[1])^2)   ; aggiornati 29/11/2021
+        err_cnt_1 = 1.0/sqrt((1.0/e1c[0])^2+(1.0/e1c[1])^2)
         ;  err_cnt_0 = stddev(p0, /nan)
         ;  err_cnt_1 = stddev(p1, /nan)
 
@@ -1346,10 +1353,13 @@ pro cal_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         cnt2Jy_0 = cal_wmean(c0c, d0c, /nan)
         cnt2Jy_1 = cal_wmean(c1c, d1c, /nan)
-        w0=1/(d0/c0)
-        w1=1/(d1/c1)
-        err_cnt2Jy_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1])
-        err_cnt2Jy_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1])
+        w0=1.0/(d0c/c0c)
+        w1=1.0/(d1c/c1c)
+        ; err_cnt2Jy_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/sqrt(w0[0]+w0[1])
+        ; err_cnt2Jy_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/sqrt(w1[0]+w1[1])
+
+        err_cnt2Jy_0 = 1.0/sqrt((1.0/d0c[0])^2+(1.0/d0c[1])^2)   ; aggiornati 29/11/2021
+        err_cnt2Jy_1 = 1.0/sqrt((1.0/d1c[0])^2+(1.0/d1c[1])^2)
 
 
         ; last check: when calibrators are not known, the resulting dummy cnt2Jy values are to be reset to -99.00

runtarget.pro

@@ -533,8 +533,10 @@ pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
           ; assigning the identified best-matching cnt2Jy values to the final attributes, for the subsequent analysis
           mylinc2j_0=linc2j_0[my_lin_interval]+mymjd*linm_0[my_lin_interval]
           mylinc2j_1=linc2j_1[my_lin_interval]+mymjd*linm_1[my_lin_interval]
-          err_mylinc2j_0=mylinc2j_0*sqrt(err_linc2j_0[my_lin_interval]^2+mymjd*err_linm_0[my_lin_interval]^2)
-          err_mylinc2j_1=mylinc2j_1*sqrt(err_linc2j_1[my_lin_interval]^2+mymjd*err_linm_1[my_lin_interval]^2)
+;          err_mylinc2j_0=mylinc2j_0*sqrt(err_linc2j_0[my_lin_interval]^2+mymjd*err_linm_0[my_lin_interval]^2) ; originale errato
+;          err_mylinc2j_1=mylinc2j_1*sqrt(err_linc2j_1[my_lin_interval]^2+mymjd*err_linm_1[my_lin_interval]^2)
+          err_mylinc2j_0=sqrt(err_linc2j_0[my_lin_interval]^2+mymjd*err_linm_0[my_lin_interval]^2)
+          err_mylinc2j_1=sqrt(err_linc2j_1[my_lin_interval]^2+mymjd*err_linm_1[my_lin_interval]^2)
         endelse
       endif
 
@@ -575,8 +577,15 @@ pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
           ; assigning the identified best-matching cnt2Jy values to the final attributes, for the subsequent analysis
           mycubc2j_0=cubc2j_0[my_cub_interval]+mymjd*cubm_0[my_cub_interval]
           mycubc2j_1=cubc2j_1[my_cub_interval]+mymjd*cubm_1[my_cub_interval]
-          err_mycubc2j_0=mycubc2j_0*sqrt(err_cubc2j_0[my_cub_interval]^2+mymjd*err_cubm_0[my_cub_interval]^2)
-          err_mycubc2j_1=mycubc2j_1*sqrt(err_cubc2j_1[my_cub_interval]^2+mymjd*err_cubm_1[my_cub_interval]^2)
+;          err_mycubc2j_0=mycubc2j_0*sqrt(err_cubc2j_0[my_cub_interval]^2+mymjd*err_cubm_0[my_cub_interval]^2)  ; originale errato
+;          err_mycubc2j_1=mycubc2j_1*sqrt(err_cubc2j_1[my_cub_interval]^2+mymjd*err_cubm_1[my_cub_interval]^2)
+          err_mycubc2j_0=sqrt(err_cubc2j_0[my_cub_interval]^2+mymjd*err_cubm_0[my_cub_interval]^2)
+          err_mycubc2j_1=sqrt(err_cubc2j_1[my_cub_interval]^2+mymjd*err_cubm_1[my_cub_interval]^2)
+          ; check per la revisione del sistema di calcolo/utilizzo dei cnt2Jy
+          ;  print, 'c2j0 medio - c2j1 medio - err orig0 - err orig1  - err c2j0 medio - err c2j1 medio
+          ;  print, mycubc2j_0, mycubc2j_1, err_cubc2j_0[my_cub_interval], err_cubc2j_1[my_cub_interval], err_mycubc2j_0, err_mycubc2j_1, format='(D6.4,1X,D6.4,1X,D6.4,1X,D6.4,1X,D6.4,1X,D6.4)'
+          ;  print, ' '
+          ; stop
         endelse
       endif
 
@@ -1122,10 +1131,13 @@ pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         peak_cnt_0 = tar_wmean(p0, e0, /nan)
         peak_cnt_1 = tar_wmean(p1, e1, /nan)
-        w0=1/(e0/p0)
-        w1=1/(e1/p1)
-        err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
-        err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
+        w0=1.0/(e0/p0)
+        w1=1.0/(e1/p1)
+        ; err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
+        ; err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
+        err_cnt_0 = 1.0/sqrt((1.0/e0[0])^2+(1.0/e0[1])^2)   ; aggiornati 29/11/2021
+        err_cnt_1 = 1.0/sqrt((1.0/e1[0])^2+(1.0/e1[1])^2)
+
         ;  err_cnt_0 = stddev(p0, /nan)
         ;  err_cnt_1 = stddev(p1, /nan)
 
@@ -1143,10 +1155,13 @@ pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
         ; XXX
         fluxdensity_0 = tar_wmean(c0, d0, /nan)
         fluxdensity_1 = tar_wmean(c1, d1, /nan)
-        w0=1/(d0/c0)
-        w1=1/(d1/c1)
-        err_fluxdensity_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1])
-        err_fluxdensity_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1])
+        w0=1.0/(d0/c0)
+        w1=1.0/(d1/c1)
+       ; err_fluxdensity_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1]) 
+       ; err_fluxdensity_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1]) 
+        err_fluxdensity_0 = 1.0/sqrt((1.0/d0[0])^2+(1.0/d0[1])^2)   ; aggiornati 29/11/2021
+        err_fluxdensity_1 = 1.0/sqrt((1.0/d1[0])^2+(1.0/d1[1])^2)
+
 
         ; last check: when calibrators are not known, the resulting dummy flux values are to be reset to -99.00
         ; as the above offset compensation affects even them (and steers them a little bit from the dummy value)
@@ -1207,10 +1222,13 @@ pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         peak_cnt_0 = tar_wmean(p0c, e0c, /nan)
         peak_cnt_1 = tar_wmean(p1c, e1c, /nan)
-        w0=1/(e0c/p0c)
-        w1=1/(e1c/p1c)
-        err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
-        err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
+        w0=1.0/(e0c/p0c)
+        w1=1.0/(e1c/p1c)
+      ;  err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1])
+      ;  err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1])
+        err_cnt_0 = 1.0/sqrt((1.0/e0c[0])^2+(1.0/e0c[1])^2)   ; aggiornati 29/11/2021
+        err_cnt_1 = 1.0/sqrt((1.0/e1c[0])^2+(1.0/e1c[1])^2)
+
         ;  err_cnt_0 = stddev(p0, /nan)
         ;  err_cnt_1 = stddev(p1, /nan)
 
@@ -1227,10 +1245,13 @@ pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, sou
 
         fluxdensity_0 = tar_wmean(c0c, d0c, /nan)
         fluxdensity_1 = tar_wmean(c1c, d1c, /nan)
-        w0=1/(d0/c0)
-        w1=1/(d1/c1)
-        err_fluxdensity_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1])
-        err_fluxdensity_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1])
+        w0=1.0/(d0c/c0c)
+        w1=1.0/(d1c/c1c)
+       ; err_fluxdensity_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1])
+       ; err_fluxdensity_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1])
+        err_fluxdensity_0 = 1.0/sqrt((1.0/d0c[0])^2+(1.0/d0c[1])^2)   ; aggiornati 29/11/2021
+        err_fluxdensity_1 = 1.0/sqrt((1.0/d1c[0])^2+(1.0/d1c[1])^2)
+
 
         ; last check: when calibrators are not known, the resulting dummy flux values are to be reset to -99.00
         ; as the above offset compensation affects even them (and steers them a little bit from the dummy value)

targetfit.pro

@@ -42,6 +42,7 @@ pro getflux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, outflux, err_o
   err_cnt=err_cnt*exp(mean(tau0)/sin(datael))
   outflux = peak_cnt*cnt2Jy
   err_outflux = outflux*sqrt((err_cnt/peak_cnt)^2+(err_cnt2Jy/cnt2Jy)^2)
+  ; print, peak_cnt, err_cnt, cnt2Jy, err_cnt2Jy
   Level=Level*cnt2Jy
   PRINTF, Unit2, " "
   PRINTF, Unit2, "outflux = (", outflux, "+-", err_outflux, ") [Jy]", format='(A,1X,D7.3,1X,A,1X,D7.3,1X,A)'