diff --git a/source/SourcePopulation.cc b/source/SourcePopulation.cc
index 7ca592596d02dd4009519ea5ade922e414bb84d9..08ebdc23422787694dfc7d986c22620d13c77e85 100644
--- a/source/SourcePopulation.cc
+++ b/source/SourcePopulation.cc
@@ -1511,7 +1511,9 @@ int SourcePopulation::analyse_sample()
                                     << lat_min1 <<" < b < "<< lat_max1 <<",  "<< lat_min2 <<" < b < "<< lat_max2<<endl;
    *txt_stream<<"converted to      region "<<long_min1a<<" < l < "<<long_max1a<<",  "<<long_min2a<<" < l < "<<long_max2a<<endl;
                                 
+   cout<<"writing summary table to results txt for "<< title<<endl;         //AWS20200622
 
+   *txt_stream<<"summary table of N(S) for "<< title <<endl;                //AWS20200622
 
    for (ii=0;ii< n_dlnN_dlnS;ii++)
    {
@@ -1536,6 +1538,41 @@ int SourcePopulation::analyse_sample()
             <<endl;
     }
 
+   // output for external plotting routines                                 //AWS20200622
+
+   *txt_stream<<"format for external plotting routines"<<endl;              //AWS20200622
+   *txt_stream<<"S = ["; 
+   for (ii=0;ii< n_dlnN_dlnS;ii++)
+       { *txt_stream<<pow(10, (lnS_min+ii*dlnS + lnS_min+(ii+1)*dlnS)/2. );
+           if(ii<n_dlnN_dlnS-1)    *txt_stream<<",";
+       }
+   *txt_stream<<"]";
+   *txt_stream<<";  idl format";
+   *txt_stream<<endl;
+
+  
+   *txt_stream<<"NS = ["; 
+   for (ii=0;ii< n_dlnN_dlnS;ii++)
+       { *txt_stream<<dlnN_dlnS[ii];
+           if(ii<n_dlnN_dlnS-1)    *txt_stream<<",";
+       }
+   *txt_stream<<"]";
+   *txt_stream<<";  idl format";
+   *txt_stream<<endl;
+
+   *txt_stream<<"set_plot,'ps'               ;idl"<<endl;
+   *txt_stream<<"device,file='"<<title<<".ps';idl"<<endl;
+   *txt_stream<<"!p.font=0                   ;idl"<<endl;
+
+   *txt_stream<<"start_plot =1 ; idl"<<endl;
+
+   // psym -6 open squares with line
+
+   *txt_stream<<"if(start_plot eq 1) then plot,S,NS,/xlog,/ylog,xrange=[1e-13,1e-5],yrange=[0.1,1e4], xstyle=1,ystyle=1, psym=-6, xtitle='S,  cm!e-2!n s!e-1!n',ytitle='N(S)', charsize=0.75,title='"<<title<<"' ; idl"<<endl;
+
+   *txt_stream<<"start_plot =0; idl";
+   *txt_stream<<"oplot,S,NS   ; idl"<<endl;
+
    *txt_stream<<"(oversampling corrected) dlnN_dlnS:  total N ="<< dlnN_dlnS_total_N  <<"           total S ="<< dlnN_dlnS_total_S<< endl;
 
    *txt_stream<<" sample region corrected for oversampling factor ="<<oversample                         <<endl;