Accessing the 6 cm radio flux density information about 
long term variability (1 year time scale) 
from the catalog (publt4.dat)
and 
information about short term variability
(1 day to 1 month) from the catalogs 
(pubcat10Jun16.dat & pubfit10Jun16)

Phil Gregory 
June 19, 2016

INTRODUCTION
     The NRAO seven-beam 14 channel receiver was used on the former Green Bank 91 m telescope, during 1986 November and 1987 October, to repeatedly survey the declination band (0° < Declination < 75°) at 4.85 GHz. This database has been used to extract variability information for sources in the GB6 catalog, a large unbiased sample of sources derived from the combined survey images. 
     This variability information consists of two parts. First, there is an indicator of the long term variability, based on two flux density measurements, one derived from the 1986 image set and the other from the 1987 image set. Long term variability information is available for 95% of the 75,162 GB6 sources.
     Second, there is an indicator of the short term variability, based on an analysis of the daily scan data that was used to construct the images. It has been possible to derive short term variability information for 20,931 GB6 sources. The exact number of daily measurements depends on source declination, ranging from a mean of 7 (both epochs) for a source at declination 0° to a mean of 20 for a declination of 73°. 
     The Green Bank 4.85 GHz survey was carried out with the following three types of data products in mind.
a) A set of sky maps ( Condon et al., AJ, 97, 1989 & AJ, 107, 1994).

b) These sky maps were used to construct catalogs of discrete radio sources ( Gregory and Condon, ApJS, 75, 1991; 
    Gregory et al.  ApJS, 103, 1996). The GB6 catalog of 75,162 discrete sources, with angular sizes < 5 arcmin and 
    flux densities > 18 mJy, was based on the sky maps contructed from both epochs. 
    
    The same software was also used to catalog 15,045 discrete sources  (-88° < Declination < -37°), with angular sizes < 15 arcmin, in the 
    Parkes-MIT-NRAO survey of the Southern Hemisphere using the same 7 feed receiver on the Parkes 64 m telescope 
    (Gregory, P.C. et al., Ap. J. Supplement 90, 173, 1994).

c) The final goal was to extract limited variability information for the entire catalog of GB6 sources. The purpose of this
    Mathematica notebook is to allow you to access this variability information. References:
    ( Gregory, Scott, and Poller, 1998, IAU Colloquium 164;  Gregory, Capak, Gasson, Scott, 2001, IAU Symposium 205 )

SETUP

Remove["Global`*"]

Off[General::"spell"]

Off[General::"spell1"]

Needs["MultivariateStatistics`"]

Needs["ComputationalGeometry`"]

Needs["ErrorBarPlots`"]

Needs["Splines`"]

Format of LongTerm Variability Catalog named publt4.dat:

publt4.dat (latest long term variability data)

1) HH:MM:SS.S        GB6 Right Ascension (B1950)
2) SS.S                     Uncertainty in right ascension (s)
3) SDD:MM:SS         GB6 Declination (B1950)
4) SS                        Uncertainty in declination (arcsec)
5) GB6 Right Ascension in decimal degrees
6) GB6 Declination in decimal degrees
7) HH:MM:SS.S        GB6 Right Ascension (J2000)
8) SDD:MM:SS         GB6 Declination (J2000)
9) GB6 Right Ascension in decimal degrees (J2000)
10) GB6 Declination in decimal degrees (J2000)
11) GB6 Galactic Latitude B in decimal degrees
12) GB6 Galactic Longitude L in decimal degrees
13) GB6 flux density (mJy)
14) GB6 flux density uncertainty (mJy)
15) Flag "E" for significant extension
16) Flag "W" for Warning
17) Flag "C"  for Confusion
18) Fitted FWHM major axis divided by point source FWHM of major axis 
19) Fitted FWHM minor axis divided by point source FWHM of minor axis
20) Fitted major axis position angle PA (degrees east of north)
21) Local sky (mJy)
22) Map pixel column number counted from left
23) Map pixel row number counted from the bottom
24) 1986 epoch flux density S86 (mJy)
25) S86 flag   ("F" = fitted flux density,   "P" = peak pixel  flux density,  "N" = no flux density)  
26) 1986 epoch flux density uncertainty  (mJy)
27) 1987 epoch flux density S87 (mJy)
28) S87 flag   ("F" = fitted flux density,   "P" = peak pixel  flux density,  "N" = no flux density) 
29) 1987 epoch flux density uncertainty  (mJy)
30) Variability index   V =  ((Abs[S86-S87])/(S86+s87))
31) Variability index error  Verr =  ((Sqrt[(((error S86))^2)+(((error S87))^2)])/(S86+S87))
32) ratio = (V/Verr)

Read in long term variability (publt4.dat)  
(you need to edit the path for the location of the catalog on your computer)

lt= Import[ "e:/Lexar/capak/longterm/publt4.dat"];

no=Length[lt]

lt1=lt[[1]]

FullForm[lt[[1]]]

ncol=Length[lt[[1]]]

Other way of inputing the data which splits the RA and DEC that contain 3 fields each separated by ":" into the 3 separate components.

(*lt = ReadList[ "e:/Lexar/capak/longterm/publt4.dat",
{Word,Word,Word,Number,Word,Word,Word,Number,Number,Number,Word,Word,Word,Word,Word,Word,Number,Number,Number,Number,Number,Number,Word,Word,Word,Number,Number,Number,Number,Number,Number,Number,Word,Number,Number,Word,Number,Number,Number,Number}, WordSeparators->{" ", ":"}];*)

Example of how to convert the String fields to Numerical fields.

(* a={"50.0","6.0",7.0} *)

(*FullForm[a]*)

(*FullForm[xx=ToExpression[a];xx]*)

(*FullForm[ToExpression[lt1]]*)

xt=Take[lt,2]

FullForm[xt]

Retrieve all Long Term (LT) source data with ratio ≥ 3 , and |B| ≤ 2.5 deg  and output source list.

The variability ratio is given by item (32) in publt4.dat which in turn depends on items (30) and (31). 
30) Variability index   V =  ((Abs[S86-S87])/(S86+s87))
31) Variability index error  Verr =  ((Sqrt[(((error S86))^2)+(((error S87))^2)])/(S86+S87))
32) ratio = (V/Verr)

lt[[1]]

lt[[1,32]]

Set lower bound on variability ratio index.

ratioMin=3;

Set upper limit on Galactic Latitude (B).

bmax=2.5;

Specify column containing Galactic Latitude (B).

bcol=11;

Below galB is a list of the B values for the selected data.

Below mj is a list of the source parameters for the selected data.

galB={};mj={};nB=0;Do[If[lt[[i,ncol]]≥ratioMin&&Abs[lt[[i,bcol]]]≤bmax,nB++;galB=Append[galB,lt[[i,bcol]]];mj=Append[mj,lt[[i]]]],{i,1,no}];nB

List of the B values for the selected data.

galB

Provide sources sorted according to increasing RA.

smg=RotateRight[Sort[RotateLeft[mj,{0,5}]],{0,5}];

(*smg=Sort[mj,OrderedQ[{#2, #1}]&];*)

nB

Choose particular source parameters to list in final table of selected sources.

printout=Table[Join[Take[mj[[i]],{7,8}],Take[mj[[i]],{11}],Take[mj[[i]],{13,14}],Take[mj[[i]],{24}],Take[mj[[i]],{26}],Take[mj[[i]],{27}],Take[mj[[i]],{29}],Take[mj[[i]],{32}]],{i,1,nB}];

Specify the header information for columns in the table of selected sources.

l1={{"RA2000","DEC2000","B","SGB6","Serr","S86","S 86err","S87","S 87err","Ratio"},{" "," "," ","mJy","mJy","mJy","mJy","mJy","mJy"," "}}

Table containing nB of the selected sources.

smj=MatrixForm[Join[l1,Take[Transpose[Transpose[printout]],nB]]]

Create the full table of the selected sources properties with header information.

data=Join[l1,printout];

To export the data file remove the (*  *) at the ends of the next command.

(*Export["e:/Lexar/capak/ratio_GalLat_selected.dat",data,"Table"];*)

Retrieve all LT variables with Flag "E" for significant extension.

lt[[1]]

FullForm[lt[[1]]]

Specify column containing normalized major axis values.

colEmajor=18;

Specify column containing normalized minor axis values.

colEminor=19;

Specify column containing significant extension flag.

colE=15;

Below colEmajor is a list of the normalized major axis values for extended sources and colEminor is a list of the normalized minor axis values for extended sources.

Below mj is a list of complete source parameters for the selected data.

majorE={};minorE={};mj={};nE=0;Do[If[lt[[i,colE]]=="E",nE++;majorE=Append[majorE,lt[[i,colEmajor]]];minorE=Append[minorE,lt[[i,colEminor]]];mj=Append[mj,lt[[i]]]],{i,1,no}];nE

List of the normalized major axes of the selected data.

majorE

Histogram of the normalized major axes of the selected data.

plEmajorHist=Histogram[majorE,PlotRange->{{1.0,3.2},All},Frame->True,ChartStyle->EdgeForm[],Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Normalized major axis",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Histogram of the normalized minor axes of the selected data.

plEminorHist=Histogram[minorE,PlotRange->{{0.5,2.5},All},Frame->True,ChartStyle->EdgeForm[],Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Normalized minor axis",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Provide sources sorted according to increasing RA.

smg=RotateRight[Sort[RotateLeft[mj,{0,5}]],{0,5}];

(*smg=Sort[mj,OrderedQ[{#2, #1}]&];*)

nE

nEsub = number of selected sources to print out.

nEsub=If[nE>50,50,nE]

Choose particular source parameters to list in final table of selected sources.

printout=Table[Join[Take[mj[[i]],{7,8}],Take[mj[[i]],{colE}],Take[mj[[i]],{13,14}],Take[mj[[i]],{colEmajor}],Take[mj[[i]],{colEminor}]],{i,1,nEsub}];

Specify the header information for columns in the table of selected sources.

l1={{"RA2000","DEC2000","Flag","(Subscript[S, GB6])","(Subscript[S, err])","Major axis","Minor axis"},{" "," "," ","mJy","mJy","",""}}

Final table of selected sources.

smj=MatrixForm[Join[l1,Take[Transpose[Transpose[printout]],nB]]]

Retrieve all LT variables with Flag "P" for peak pixel flux density 
in 86 or 87 instead of fitted flux density.

lt[[1]]

FullForm[lt[[1]]]

lt[[1,29]]

Specify columns containing S86  value and error.

colS86=24;colSerr86=26;

Specify columns containing S87  value and error.

colS87=27;colSerr87=29;

Specify columns containing S86 "P" flag &  S87 "P" flag.

colP86=25;colP87=28;

Below mj is a list of complete source parameters for the selected data.

mj={};nP=0;nP86=0;nP87=0;Do[If[lt[[i,colP86]]=="P",nP86++];If[lt[[i,colP87]]=="P",nP87++];If[lt[[i,colP86]]=="P"||lt[[i,colP87]]=="P",nP++;mj=Append[mj,lt[[i]]]],{i,1,no}];nP
nP86
nP87

n86noS = Number of 86 sources for which no flux density estimate possible
                 (ie., S86 = -999).

s86=Transpose[mj][[colS86]];n86noS=0;Do[If[s86[[i]]==-999,n86noS++],{i,1,nP}];n86noS

n87noS = Number of 87 sources for which no flux density estimate possible
                 (ie., S87 = -999).

s87=Transpose[mj][[colS87]];n87noS=0;Do[If[s87[[i]]==-999,n87noS++],{i,1,nP}];n87noS

Histogram of GB 6 flux densities for the selected data.

gb6S=Transpose[mj][[13]];

plgb6SHist=Histogram[gb6S,PlotRange->{{0,90},All},Frame->True,ChartStyle->EdgeForm[],Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["GB6 flux density",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Provide sources sorted according to increasing RA.

smg=RotateRight[Sort[RotateLeft[mj,{0,5}]],{0,5}];

(*smg=Sort[mj,OrderedQ[{#2, #1}]&];*)

nP

nPsub = number of selected sources to print out.

nPsub=If[nP>50,50,nP]

Choose particular source parameters to list in final table of selected sources.

printout=Table[Join[Take[mj[[i]],{7,8}],Take[mj[[i]],{colP86}],Take[mj[[i]],{colP87}],Take[mj[[i]],{13,14}],Take[mj[[i]],{colS86}],Take[mj[[i]],{colS87}]],{i,1,nPsub}];

Specify the header information for columns in the table of selected sources.

l1={{"RA2000","DEC2000","P Flag86","P Flag87","(Subscript[S, GB6])","(Subscript[S, err])","S86","S87"},{" "," "," "," ","mJy","mJy","mJy","mJy"}}

Final table of selected sources.

smj=MatrixForm[Join[l1,Take[Transpose[Transpose[printout]],nB]]]

Retrieve all LT variables within θ arcmin of a J2000 reference coordinate.
Uses command StringSplit  (faster than StringToStream)

Specify θ size of acceptance window in arcmin

θ=30.0

θrad=θ*(π/(60*180))

Specify reference coordinate RA and DEC.

raRef="09:36:33.1"

decRef="+36:49:10"

Convert raRef in "HH:MM:SS.S" to hms, a list of three Numbers {HH,MM,SS.S}

hms=ToExpression[StringSplit[raRef,":"]]

FullForm[hms]

Convert hms to radians

raRefRad=(hms[[1]]*15.0+hms[[2]]/4.0+hms[[3]]/240.0)*(π/180)

Convert raRef in "HH:MM:SS.S" to a list of three Numbers {HH,MM,SS.S}

dms=ToExpression[StringSplit[decRef,":"]]

Convert dms to radians

decRefRad=(dms[[1]]+dms[[2]]/60.0+dms[[3]]/3600.0)*(π/180)

no

Specifies the columns that contain J2000 RA and DEC.

raJ2000Col=7;decJ2000Col=8;

Portion of code that finds all sources (output as lis) that are within an angular distance of θrad from the reference position. The angular distance of the select sources are output as anglis.

timerun=AbsoluteTiming[
lis={};anglis={};
Do[raX=lt[[i,raJ2000Col]];decX=lt[[i,decJ2000Col]];hmsX=ToExpression[StringSplit[raX,":"]];raXRad=(hmsX[[1]]*15.0+hmsX[[2]]/4.0+hmsX[[3]]/240.0)*(π/180);
dmsX=ToExpression[StringSplit[decX,":"]];decXRad=(dmsX[[1]]+dmsX[[2]]/60.0+dmsX[[3]]/3600.0)*(π/180);
ang=ArcCos[Sin[decXRad]*Sin[decRefRad]+Cos[decXRad]*Cos[decRefRad]*Cos[(raXRad-raRefRad)]];If[ang⩽θrad,lis=Append[lis,lt[[i]]];anglis=Append[anglis,ang*180*60/π]],{i,1,no}
];nsel=Length[lis]
]

Separations of selected sources in arcmin from reference position.

anglis

Selected source parameters for reference position.

lis

nPsub = number of selected sources to print out.

nAsub=If[nsel>50,50,nsel]

Choose particular source parameters to list in final table of selected sources.

Specify columns containing S86  value and error.

colS86=24;colSerr86=26;

Specify columns containing S87  value and error.

colS87=27;colSerr87=29;

printout=Table[Join[Take[lis[[i]],{7,8}],{anglis[[i]]},Take[lis[[i]],{13,14}],Take[lis[[i]],{colS86}],Take[lis[[i]],{colSerr86}],Take[lis[[i]],{colS87}],Take[lis[[i]],{colSerr87}]],{i,1,nAsub}];

Specify the header information for columns in the table of selected sources.

l1={{"RA2000","DEC2000","Angle","(Subscript[S, GB6])","(Subscript[S, err])","S86","(Subscript[S, 86  err])","S87","(Subscript[S87, err])"},{" "," ","Arc min ","mJy","mJy","mJy","mJy","mJy","mJy"}}

Final table of selected sources.

smj=MatrixForm[Join[l1,Take[Transpose[Transpose[printout]],nAsub]]]

Plot a Parabolic Projection In Celestial Coordinates (J2000) with ratio ≥3.

RA defined positive measured eastward from the vernal equinox

k=(π/180);

xc[ϕc_,θc_]:=-ϕc (2 Cos[((2 k θc)/3)] - 1)

yc[ϕc_,θc_]:=180 Sin[((k θc)/3)]

N[xc[180,0]]

l1=ParametricPlot[{xc[-180,θc],yc[-180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l2=ParametricPlot[{xc[-120,θc],yc[-120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l3=ParametricPlot[{xc[-60,θc],yc[-60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l4=ParametricPlot[{xc[60,θc],yc[60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l5=ParametricPlot[{xc[120,θc],yc[120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l6=ParametricPlot[{xc[180,θc],yc[180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];

b1=ParametricPlot[{xc[ϕc,60],yc[ϕc,60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b2=ParametricPlot[{xc[ϕc,30],yc[ϕc,30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b3=ParametricPlot[{xc[ϕc,-30],yc[ϕc,-30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b4=ParametricPlot[{xc[ϕc,-60],yc[ϕc,-60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];

Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4}]

First entry in the long term variability data set (publt4.dat)

lt[[1]]

The variability ratio is given by item (32) in publt4.dat which in turn depends on items (30) and (31). 
30) Variability index   V =  ((Abs[S86-S87])/(S86+s87))
31) Variability index error  Verr =  ((Sqrt[(((error S86))^2)+(((error S87))^2)])/(S86+S87))
32) ratio = (V/Verr)

lt[[1]][[32]]

Select sources with ratio >3.0

nsel=0;vg={};Do[If[lt[[i]][[32]]>3.0,nsel++;vg=Append[vg,{lt[[i]][[9]],lt[[i]][[10]]}],],{i,1,no}];nsel

Construct a table of RA, DEC pairs for each selected source

y=Table[{xc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]],yc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]]},{i,1,nsel}];

Plot the locations of the selected sources in RA, DEC

p1=ListPlot[y,PlotRange->{{-180,180},{-90,90}}];

radecpl=Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4,p1},DisplayFunction->$DisplayFunction]

Format of 2 Short Term Variability Catalog files:
 Intended to provide coarse information on short term variability within each of the 1986 and 1987 
 observing epochs. Based on a Bayesian analysis of the 7 feed data for each day where the beams are  separated by 1 HPBW in the direction perpendicular to the sky track.

pubcat10Jun16.dat (short term summary data with variability statistics for daily flux density detections)

Format of  "pubcat10Jun16.dat" that contains summary information and statistics on the daily flux density measurements.
1)   Serial Number    014050196
2)    Δθ = (daily scan data derived position - GB6 position)  arcmin                     
3)   HH:MM:SS.S        GB6 Right Ascension (RA) (B1950)
4)   SS.S                     Uncertainty in right ascension (s)
5)   SDD:MM:SS         GB6 Declination (Dec) (B1950)
6)   SS                        Uncertainty in declination (arcsec)
7)    GB6 RA      (B1950, decimal degrees)
8)    GB6 Dec    (B1950, decimal degrees)
9)     HH:MM:SS.S        GB6 RA   (J2000) 
10)   SDD:MM:SS.SS  GB6 Dec  (J2000) 
11)  GB6 RA     (J2000, decimal degrees)
12)  GB6 Dec   (J2000, decimal degrees)        
13)  GB6 Flux density            (mJy)
14)  GB6 Flux density Error    (mJy)
15)  Number of confusing sources subtracted
16)  Maximum Scan Flux density    (mJy)
17)  Fitted RA in minutes of time        (1950)
18)  Fitted RA in minutes of time error 
19)  Fitted Dec in decimal degrees    (1950)
20)  Fitted Dec in decimal degrees error  
21)  Galactic Latitude       (decimal degrees)
22)  Galactic Longitude    (decimal degrees)
23)  Number of fits in epoc 1986
24)  Number of fits in epoc 1987
25)  Average Flux density (both epochs)      (mJy)
26)  Average Flux density epoc 1986            (mJy)
27)  Average Flux density epoc 1987            (mJy)
28)  Weighted average Flux density (both epochs)    (mJy)
29)    Error in weighted average                                  (mJy)
30)  Weighted average Flux density (OverBar[(Subscript[S, wt86])]) epoc 1986          (mJy)
31)    Error in weighted average1986                          (mJy)
32)  Weighted average Flux density (OverBar[(Subscript[S, wt87])]) epoc 1987          (mJy)
33)    Error in weighted average 1987                         (mJy)
34)  Modulation Index  epoc 1986  given by (1/(OverBar[(Subscript[S, wt86])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i86])-(OverBar[(Subscript[S, wt86])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))
35)  Modulation Index  epoc 1987 given by (1/(OverBar[(Subscript[S, wt87])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i87])-(OverBar[(Subscript[S, wt87])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))  
36)  Modulation Index due to noise from 1000 simulations for epoc 1986   (-9 means ignore)
37)  Modulation Index due to noise from 1000 simulations for epoc 1987   (-9 means ignore)
38)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1986   (-9 means ignore)
39)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1987   (-9 means ignore) 
40)  Corrected Modulation Index epoc 1986
41)  Corrected Modulation Index epoc 1987
42)  -(Subscript[log, 10]) (p-value)  
43)  -(Subscript[log, 10]) (p-value)  for epoc  1   
44)  -(Subscript[log, 10]) (p-value) for epoc 2

Explanation of the - Log(p-value) variability indicator. We employ the common  (χ^2) statistic in a frequentist statistical test of the null hypothesis that the source is a constant, equal to the weighted mean of the daily flux densities. The p-value measures the area in the tail of a (χ^2) distribution with m degrees of freedom for a value of (χ^2) greater than or equal to the measured (χ^2) for the observed set of daily flux densities. The number of degrees of freedom m = N - 1, where N is the number of daily flux densities. The total area under the (χ^2) distribution is 1.

Suppose the area in the tail region is 0.05 or 5%. What does this mean? It means that if the flux density of this object is really constant and we repeately obtained sets of N measurements under the same conditions (hypothetical data sets), then only 5 % of the (χ^2) values derived from these hypothetical data sets would be expected to be greater than or equal to the value of (χ^2) computed for the one actually measured data set. The statistical confidence for rejecting the null hypothesis that the source has a constant flux density is = 1 - (p-value), which in this example is 0.95 or 95 %. Thus, the smaller the p-value the more likely the source is not constant. It is convenient to compute the statistic -(Subscript[log, 10]) (p-value). The larger the value of -(Subscript[log, 10]) (p-value) the more confident we are in rejecting the null hypothesis that the source has a constant flux density.

pubfit10Jun16.dat (daily flux density data)

Format of  "pubfit10Jun16.dat" data. 
1) Unique source serial number
2)  Scan number                     
3)  Daily flux density  (mJy)
4)  Daily flux density uncertainty  (mJy)        
5)  Modified Julian Day number (accurate to 4 decimal places)   MJD = JD - 2,400,000.5 
6)  Flag   (0= fit but (ave or wtd ave) flux too low compared to GB6, 1= successful fit)

Read in Short Term (ST) variability summary catalog with statistics (pubcat10Jun16.dat)
(you need to edit the path for the location of the catalog on your computer)

Format of  "pubcat10Jun16.dat" and earlier version "pubcat.dat".
1)   Serial Number    014050196
2)    Δθ = (daily scan data derived position - GB6 position)  arcmin                     
3)   HH:MM:SS.S        GB6 Right Ascension (RA) (B1950)
4)   SS.S                     Uncertainty in right ascension (s)
5)   SDD:MM:SS         GB6 Declination (Dec) (B1950)
6)   SS                        Uncertainty in declination (arcsec)
7)    GB6 RA      (B1950, decimal degrees)
8)    GB6 Dec    (B1950, decimal degrees)
9)     HH:MM:SS.S        GB6 RA   (J2000) 
10)   SDD:MM:SS.SS  GB6 Dec  (J2000) 
11)  GB6 RA     (J2000, decimal degrees)
12)  GB6 Dec   (J2000, decimal degrees)        
13)  GB6 Flux density            (mJy)
14)  GB6 Flux density Error    (mJy)
15)  Number of confusing sources subtracted
16)  Maximum Scan Flux density    (mJy)
17)  Fitted RA in minutes of time        (1950)
18)  Fitted RA in minutes of time error 
19)  Fitted Dec in decimal degrees    (1950)
20)  Fitted Dec in decimal degrees error  
21)  Galactic Latitude       (decimal degrees)
22)  Galactic Longitude    (decimal degrees)
23)  Number of fits in epoc 1986
24)  Number of fits in epoc 1987
25)  Average Flux density (both epochs)      (mJy)
26)  Average Flux density epoc 1986            (mJy)
27)  Average Flux density epoc 1987            (mJy)
28)  Weighted average Flux density (both epochs)    (mJy)
29)    Error in weighted average                                  (mJy)
30)  Weighted average Flux density epoc 1986          (mJy)
31)    Error in weighted average1986                          (mJy)
32)  Weighted average Flux density epoc 1987          (mJy)
33)    Error in weighted average 1987                         (mJy)
34)  Modulation Index  epoc 1986  given by (1/(OverBar[(Subscript[S, wt86])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i86])-(OverBar[(Subscript[S, wt86])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))
35)  Modulation Index  epoc 1987 given by (1/(OverBar[(Subscript[S, wt87])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i87])-(OverBar[(Subscript[S, wt87])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))  
36)  Modulation Index due to noise from 1000 simulations for epoc 1986   (-9 means ignore)
37)  Modulation Index due to noise from 1000 simulations for epoc 1987   (-9 means ignore)
38)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1986   (-9 means ignore)
39)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1987   (-9 means ignore) 
40)  Corrected Modulation Index epoc 1986
41)  Corrected Modulation Index epoc 1987
42)  -(Subscript[log, 10]) (p-value)  
43)  -(Subscript[log, 10]) (p-value)  for epoc  1   
44)  -(Subscript[log, 10]) (p-value) for epoc 2

st= Import[ "e:/Lexar/capak/pubcat10Jun16.dat"];

ns=Length[st]

FullForm[st[[1]]]

nscol=Length[st[[1]]]

Examine maximum Scan Flux (mJy).

Min[Transpose[st][[16]]]

Max[Transpose[st][[16]]]

Median[Transpose[st][[16]]]

(*Export["e:/Lexar/capak/longterm/xt1_temp.dat",xt,"Table",PageWidth->Infinity];*)

Compare weighted average of pubcat10Jun16.dat daily flux densities (both epochs) to GB 6 flux densities.

Read in constsrc.dat (unresolved steep spectral index constant flux density sources).

List of constant source ID's.

idc={2040339,1404039,1404042,1404070,14040121,2204044,3004038,3004040,30040108,30040146,30040202,7404075,74040184,74040225,74040258,74040275,74040283,74040303,74040325,74040407,74040412,9404048,9404053,9404054,9404061,94040121,94040233,94040266,94040268,94040295,94040305,94040319,94040336,11404043,11404044,11404047,11404097,114040105,114040127,114040143,114040145,114040156,114040171,114040172,114040184,12204056,12204065,12204073,12204096,122040109,122040190,122040251,122040273,14204026,14204057,142040108,142040115,142040239,142040285,142040289,142040293,162040156,162040178,162040208,162040210,162040217,162040323,162040355,162040372,162040395,17004043,162040421,17004096,170040109,14020154,50020150,6205078,70050124,82050125,90050310,94010242,94010267,10201080,102060204,114070105,150010327,17000046,170060159,194060155,230000138};

noc=Length[idc]

The first source information in pubcat--.dat data.

First[st]

Extract the pubcat.dat data  corresponding to the constant source ID's.

AbsoluteTiming[cst={};
Do[id=idc[[i]];iCon=Select[st,#[[1]]==id&];If[Length[iCon]>0,cst=Append[cst,iCon[[1]]]],{i,1,noc}
];ncst=Length[cst]]

Last[cst]

Extract the GB6 flux densities and corresponding weighted daily averages.

tcst=Transpose[cst];sComp=Transpose[{tcst[[13]],tcst[[28]]}];Length[sComp]

Uncertainties in weighted daily averages.

derr=tcst[[29]]

Sorted pairs of GB6 flux densities and weighted daily averages.

sxy=Sort[sComp]

(*dataerr=Table[{sComp[[i]],ErrorBar[tcsp[[29]][[i]]]},{i,noc}];*)

plsCon=ListLogLogPlot[sComp,Frame->True,PlotRange->All,Axes->False,AspectRatio->0.5,PlotStyle->Black,PlotLabel->Style["",{"Helvetica-Bold",11}],FrameLabel->{Style["(Subscript[Log, 10]) GB6  S (mJy)",{"Helvetica-Bold",14}],
		Style["(Subscript[Log, 10]) Weighted daily ave.  S (mJy)",{"Helvetica-Bold",14}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Compute joint posterior density for slope η and intercept α for robust linear fit of daily average flux

tsxy=Transpose[sxy];

tsxy[[1]]

Note in this analysis I would have liked to have based the error term on sqrt of the sum of the squares of the error in wtd ave flux density and GB6 flux density but unfortunately I didn't keep enough significant figures in the sum of the weights term. However,
I tried various extreme assumptions like a fixed error size to see what effect it had. Below we use the Absolute magnitude normalized difference and following that the usual Gaussian noise assumption.

fn1=Exp[-Plus@@Abs[(((tcst[[28]]-α-η*tcst[[13]]))/(Sqrt[((tcst[[14]])^2)+((tcst[[29]])^2)]))]]/Exp[-Plus@@Abs[(((tcst[[28]]-0.0-1.0*tcst[[13]]))/(Sqrt[((tcst[[14]])^2)+((tcst[[29]])^2)]))]];

ContourPlot[fn1,{α,-6,6},{η,.80,1.2},PlotRange->All,PlotPoints->60,ContourShading->False,PlotLabel->Style["(B)  Probability Contours",{"Helvetica-Bold",12}],FrameLabel->{Style["α axis",{"Helvetica-Bold",12}],
		Style["η axis",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Best linear fit is not significantly different form a slope of 1.0 .

sd1=0.0+1.0*sgb6;

le=LogLogPlot[sd1,{sgb6,20,12000},PlotRange->All,PlotStyle->Black]

Show[plsCon,le]

Plot weighted daily average flux with error versus GB6 flux

yw1=Table[{Log[10,cst[[i]][[13]]],Log[10,cst[[i]][[28]]],Log[10,(cst[[i]][[28]]-cst[[i]][[29]])],Log[10,(cst[[i]][[28]]+cst[[i]][[29]])]},{i,1,noc}];

tary=Transpose[yw1];

xypt=Transpose[{N[tary[[1]]],N[tary[[2]]]}];

xpt=N[tary[[1]]];

ler=N[tary[[3]]];

uer=N[tary[[4]]];

xyplt=ListPlot[xypt,PlotRange->{{1.2,4.3},{1.2,4.3}},PlotStyle->Black,Frame->True,FrameLabel->{Style["(Subscript[Log, 10]) GB6  S (mJy)",{"Helvetica-Bold",14}],
		Style["(Subscript[Log, 10]) Weighted daily ave.  S (mJy)",{"Helvetica-Bold",14}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Plot lines between pairs of ler and uer for use as error bars .

errplt=Show[Graphics[{Thickness[0.003],Table[{Line[{{xpt[[i]],ler[[i]]},{xpt[[i]],uer[[i]]}}]},
				{i,1,noc}]}],PlotRange->{{1.2,4.3},{1.2,4.3}},Frame->True]

Create a st. line of slope 1.0 for overplotting.

lin={{Log[10,10.],Log[10,10.0]},{Log[10,18000.],Log[10,18000.]}}

le=ListPlot[lin,PlotStyle->{Black},PlotRange->All,Joined->True]

Plot points, errbars and st. line.

flx=Show[xyplt,errplt,le]

ph=Histogram[((N[tcst[[28]]]-N[tcst[[13]]])/(Sqrt[((N[tcst[[14]]])^2)+((N[tcst[[29]]])^2)])),{-4.5,4.5,0.5},PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(wtd\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

pg=Plot[noc*0.5*PDF[NormalDistribution[0,1],x],{x,-4.5,4.5},Frame->True]

plHistDeltaS=Show[ph,pg,PlotRange->{{-4.5,4.5},All},Frame->True,FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(wtd\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",12}],
		Style["Number",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Compare weighted average of pubcat10Jun16.dat daily flux densities (1986 epoch) to GB 6 flux densities.

Read in constsrc.dat (unresolved steep spectral index constant flux density sources).

List of constant source ID's.

idc={2040339,1404039,1404042,1404070,14040121,2204044,3004038,3004040,30040108,30040146,30040202,7404075,74040184,74040225,74040258,74040275,74040283,74040303,74040325,74040407,74040412,9404048,9404053,9404054,9404061,94040121,94040233,94040266,94040268,94040295,94040305,94040319,94040336,11404043,11404044,11404047,11404097,114040105,114040127,114040143,114040145,114040156,114040171,114040172,114040184,12204056,12204065,12204073,12204096,122040109,122040190,122040251,122040273,14204026,14204057,142040108,142040115,142040239,142040285,142040289,142040293,162040156,162040178,162040208,162040210,162040217,162040323,162040355,162040372,162040395,17004043,162040421,17004096,170040109,14020154,50020150,6205078,70050124,82050125,90050310,94010242,94010267,10201080,102060204,114070105,150010327,17000046,170060159,194060155,230000138};

noc=Length[idc]

First entry in pubcat--.dat.

First[st]

Extract the pubcat2.dat data  corresponding to the constant source ID's.

AbsoluteTiming[cst={};
Do[id=idc[[i]];iCon=Select[st,#[[1]]==id&];If[Length[iCon]>0,cst=Append[cst,iCon[[1]]]],{i,1,noc}
];ncst=Length[cst]]

Last[cst]

Extract the GB6 flux densities and corresponding weighted daily averages.

tcst=Transpose[cst];sCompErr=Transpose[{tcst[[13]],tcst[[14]],tcst[[30]],tcst[[31]]}];Length[sCompErr]

sCompErr

sCompErrGood=Select[sCompErr,#[[4]]>1.0&]

Length[sCompErrGood]

tsG=Transpose[sCompErrGood];

sComp=Transpose[{tsG[[1]],tsG[[3]]}]

n86=Length[sComp]

Uncertainties in weighted daily averages.

derr=tsG[[4]]

Sorted pairs of GB6 flux densities and weighted daily averages.

sxy=Sort[sComp]

(*dataerr=Table[{sComp[[i]],ErrorBar[tcsp[[29]][[i]]]},{i,noc}];*)

plsCon=ListLogLogPlot[sComp,Frame->True,PlotRange->All,Axes->False,AspectRatio->0.5,PlotStyle->Black,PlotLabel->Style["",{"Helvetica-Bold",11}],FrameLabel->{Style["(Subscript[Log, 10]) GB6  S (mJy)",{"Helvetica-Bold",14}],
		Style["(Subscript[Log, 10]) Weighted daily ave. (Subscript[S, 86]) (mJy)",{"Helvetica-Bold",14}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Compare with a slope of 1.0 and intercept of 0.0 .

sd1=0.0+1.0*sgb6;

le=LogLogPlot[sd1,{sgb6,20,12000},PlotRange->All,PlotStyle->Black]

Show[plsCon,le]

Plot weighted 1986 daily average flux with error versus GB6 flux

tsG[[1]]

yw1=Table[{Log[10,tsG[[1]][[i]]],Log[10,tsG[[3]][[i]]],Log[10,(tsG[[3]][[i]]-tsG[[4]][[i]])],Log[10,(tsG[[3]][[i]]+tsG[[4]][[i]])]},{i,1,n86}];

tary=Transpose[yw1];

xypt=Transpose[{N[tary[[1]]],N[tary[[2]]]}];

xpt=N[tary[[1]]];

ler=N[tary[[3]]];

uer=N[tary[[4]]];

xyplt=ListPlot[xypt,PlotRange->{{1.2,4.3},{1.2,4.3}},PlotStyle->Black,Frame->True,FrameLabel->{Style["(Subscript[Log, 10]) GB6  S (mJy)",{"Helvetica-Bold",14}],
		Style["(Subscript[Log, 10]) Weighted daily ave. (Subscript[S, 86]) (mJy)",{"Helvetica-Bold",14}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Plot lines between pairs of ler and uer for use as error bars .

errplt=Show[Graphics[{Thickness[0.003],Table[{Line[{{xpt[[i]],ler[[i]]},{xpt[[i]],uer[[i]]}}]},
				{i,1,n86}]}],PlotRange->{{1.2,4.3},{1.2,4.3}},Frame->True]

Create a st. line of slope 1.0 for overplotting.

lin={{Log[10,10.],Log[10,10.0]},{Log[10,18000.],Log[10,18000.]}}

le=ListPlot[lin,PlotStyle->{Black},PlotRange->All,Joined->True]

Plot points, errbars and st. line.

flx=Show[xyplt,errplt,le]

ph=Histogram[((N[tcst[[30]]]-N[tcst[[13]]])/(Sqrt[((N[tcst[[14]]])^2)+((N[tcst[[31]]])^2)])),{-4.5,4.5,0.5},PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(wtd\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

pg=Plot[n86*0.5*PDF[NormalDistribution[0,1],x],{x,-4.5,4.5},Frame->True]

plHistDeltaS=Show[ph,pg,PlotRange->{{-4.5,4.5},All},Frame->True,FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(wtd\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",12}],
		Style["Number",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Compare weighted average of pubcat10Jun16.dat daily flux densities (1987 epoch) to GB 6 flux densities.

Read in constsrc.dat (unresolved steep spectral index constant flux density sources).

List of constant source ID's.

idc={2040339,1404039,1404042,1404070,14040121,2204044,3004038,3004040,30040108,30040146,30040202,7404075,74040184,74040225,74040258,74040275,74040283,74040303,74040325,74040407,74040412,9404048,9404053,9404054,9404061,94040121,94040233,94040266,94040268,94040295,94040305,94040319,94040336,11404043,11404044,11404047,11404097,114040105,114040127,114040143,114040145,114040156,114040171,114040172,114040184,12204056,12204065,12204073,12204096,122040109,122040190,122040251,122040273,14204026,14204057,142040108,142040115,142040239,142040285,142040289,142040293,162040156,162040178,162040208,162040210,162040217,162040323,162040355,162040372,162040395,17004043,162040421,17004096,170040109,14020154,50020150,6205078,70050124,82050125,90050310,94010242,94010267,10201080,102060204,114070105,150010327,17000046,170060159,194060155,230000138};

noc=Length[idc]

Extract the pubcat2.dat data  corresponding to the constant source ID's.

First[st]

AbsoluteTiming[cst={};
Do[id=idc[[i]];iCon=Select[st,#[[1]]==id&];If[Length[iCon]>0,cst=Append[cst,iCon[[1]]]],{i,1,noc}
];ncst=Length[cst]]

Last[cst]

Extract the GB6 flux densities and corresponding weighted daily averages.

tcst=Transpose[cst];sCompErr=Transpose[{tcst[[13]],tcst[[14]],tcst[[32]],tcst[[33]]}];Length[sCompErr]

sCompErrGood=Select[sCompErr,#[[4]]>0&]

tsG=Transpose[sCompErrGood];

sComp=Transpose[{tsG[[1]],tsG[[3]]}]

n87=Length[sComp]

Uncertainties in weighted daily averages.

derr=tsG[[4]]

Sorted pairs of GB6 flux densities and weighted daily averages.

sxy=Sort[sComp]

(*dataerr=Table[{sComp[[i]],ErrorBar[tcsp[[29]][[i]]]},{i,noc}];*)

plsCon=ListLogLogPlot[sComp,Frame->True,PlotRange->All,Axes->False,AspectRatio->0.5,PlotStyle->Black,PlotLabel->Style["",{"Helvetica-Bold",11}],FrameLabel->{Style["(Subscript[Log, 10]) GB6  S (mJy)",{"Helvetica-Bold",14}],
		Style["(Subscript[Log, 10]) Weighted daily ave. (Subscript[S, 87]) (mJy)",{"Helvetica-Bold",14}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Compare with a slope of 1.0 and intercept of 0.0 .

sd1=0.0+1.0*sgb6;

le=LogLogPlot[sd1,{sgb6,20,12000},PlotRange->All,PlotStyle->Black]

Show[plsCon,le]

Plot weighted 1987 daily average flux with error versus GB6 flux

tsG[[1]]

yw1=Table[{Log[10,tsG[[1]][[i]]],Log[10,tsG[[3]][[i]]],Log[10,(tsG[[3]][[i]]-tsG[[4]][[i]])],Log[10,(tsG[[3]][[i]]+tsG[[4]][[i]])]},{i,1,n87}];

tary=Transpose[yw1];

xypt=Transpose[{N[tary[[1]]],N[tary[[2]]]}];

xpt=N[tary[[1]]];

ler=N[tary[[3]]];

uer=N[tary[[4]]];

xyplt=ListPlot[xypt,PlotRange->{{1.2,4.3},{1.2,4.3}},PlotStyle->Black,Frame->True,FrameLabel->{Style["(Subscript[Log, 10]) GB6  S (mJy)",{"Helvetica-Bold",14}],
		Style["(Subscript[Log, 10]) Weighted daily ave. (Subscript[S, 87]) (mJy)",{"Helvetica-Bold",14}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Plot lines between pairs of ler and uer for use as error bars .

errplt=Show[Graphics[{Thickness[0.003],Table[{Line[{{xpt[[i]],ler[[i]]},{xpt[[i]],uer[[i]]}}]},
				{i,1,n87}]}],PlotRange->{{1.2,4.3},{1.2,4.3}},Frame->True]

Create a st. line of slope 1.0 for overplotting.

lin={{Log[10,10.],Log[10,10.0]},{Log[10,18000.],Log[10,18000.]}}

le=ListPlot[lin,PlotStyle->{Black},PlotRange->All,Joined->True]

Plot points, errbars and st. line.

flx=Show[xyplt,errplt,le]

ph=Histogram[((N[tcst[[32]]]-N[tcst[[13]]])/(Sqrt[((N[tcst[[14]]])^2)+((N[tcst[[33]]])^2)])),{-4.5,4.5,0.5},PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(wtd\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

pg=Plot[n87*0.5*PDF[NormalDistribution[0,1],x],{x,-4.5,4.5},Frame->True]

plHistDeltaS=Show[ph,pg,PlotRange->{{-4.5,4.5},All},Frame->True,FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(wtd\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",12}],
		Style["Number",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Histogram ((((Subscript[S, daily]) - (Subscript[S, GB6])))/(Subscript[S, err])) for constant sources (both epochs)

Read list of constant source ID's (unresolved steep spectral index constant flux density sources).

idc={2040339,1404039,1404042,1404070,14040121,2204044,3004038,3004040,30040108,30040146,30040202,7404075,74040184,74040225,74040258,74040275,74040283,74040303,74040325,74040407,74040412,9404048,9404053,9404054,9404061,94040121,94040233,94040266,94040268,94040295,94040305,94040319,94040336,11404043,11404044,11404047,11404097,114040105,114040127,114040143,114040145,114040156,114040171,114040172,114040184,12204056,12204065,12204073,12204096,122040109,122040190,122040251,122040273,14204026,14204057,142040108,142040115,142040239,142040285,142040289,142040293,162040156,162040178,162040208,162040210,162040217,162040323,162040355,162040372,162040395,17004043,162040421,17004096,170040109,14020154,50020150,6205078,70050124,82050125,90050310,94010242,94010267,10201080,102060204,114070105,150010327,17000046,170060159,194060155,230000138};

noc=Length[idc]

Extract the pubcat--.dat data  corresponding to the constant source ID's.

AbsoluteTiming[cst={};
Do[id=idc[[i]];iCon=Select[st,#[[1]]==id&];If[Length[iCon]>0,cst=Append[cst,iCon[[1]]]],{i,1,noc}
];ncst=Length[cst]]

Last[cst]

Here are the GB6 flux densities.

tcst=Transpose[cst];tcst[[13]]

Extract the corresponding daily flux densities and errors.

Read in Short Term (ST) daily flux values (pubfit10Jun16.dat)

dt=Import[ "e:/Lexar/capak/pubfit10Jun16.dat"];Length[dt]

Format of  "pubfit10Jun16.dat" data. 
1) Unique source serial number
2)  Scan number                     
3)  Daily flux density  (mJy)
4)  Daily flux density uncertainty  (mJy)        
5)  Modified Julian Day number (accurate to 4 decimal places)   MJD = JD - 2,400,000.5 
6)  Flag   (0= fit but (ave or wtd ave) flux too low compared to GB6, 1= successful fit)

Length[dt[[1]]]

dt[[1]]

dti=Select[dt,#[[1]]==idc[[20]]&]

ndti=Length[dti]

tdti=Transpose[dti]

tdti[[3]]

tdti[[4]]

Here is the pubcat.dat data for source 20 with the GB 6 flux density and error in columns 13 & 14, respectively.

cst[[20]]

Here is the GB6 flux density.

tcst[[13]][[20]]

Here is the GB6 flux density error.

tcst[[14]][[20]]

dSE=(tdti[[3]]-tcst[[13]][[20]])

(Sqrt[((tdti[[4]])^2)+(((tcst[[14]][[1]]))^2)])

(tdti[[3]]-tcst[[13]][[20]])/(Sqrt[((tdti[[4]])^2)+((tcst[[14]][[20]])^2)])

lisDF={};Do[dti=Select[dt,#[[1]]==idc[[i]]&];ndti=Length[dti];If[ndti>0,tdti=Transpose[dti];dSE=(tdti[[3]]-tcst[[13]][[i]])/(Sqrt[((tdti[[4]])^2)+((tcst[[14]][[i]])^2)]);lisDF=Append[lisDF,dSE]],{i,1,noc}];Length[lisDF]

nlisDF=Length[Flatten[lisDF]]

dh=0.25;

ph=Histogram[Flatten[lisDF],{-5,5,dh},PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,ChartStyle->Gray,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(daily\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

N[Integrate[PDF[NormalDistribution[0,1],x],{x,-5,5}]]

pg=Plot[nlisDF*dh*PDF[NormalDistribution[0,1],x],{x,-5,5},PlotStyle->{Black,Thick},Frame->True]

plHistDeltaS=Show[ph,pg,PlotRange->{{-5,5},All},Frame->True,FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(daily\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",12}],
		Style["Number",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

(*Export["e:/Lexar/capak/examineVarLT_ST_19Jun16M7_HistDailyDeltaS_To_Error.eps",plHistDeltaS,ImageSize->72*7];*)

Histogram ((((Subscript[S, daily]) - (Subscript[S, GB6])))/(Subscript[S, err])) for constant sources (epoch 1986)

Read list of constant source ID's (unresolved steep spectral index constant flux density sources).

idc={2040339,1404039,1404042,1404070,14040121,2204044,3004038,3004040,30040108,30040146,30040202,7404075,74040184,74040225,74040258,74040275,74040283,74040303,74040325,74040407,74040412,9404048,9404053,9404054,9404061,94040121,94040233,94040266,94040268,94040295,94040305,94040319,94040336,11404043,11404044,11404047,11404097,114040105,114040127,114040143,114040145,114040156,114040171,114040172,114040184,12204056,12204065,12204073,12204096,122040109,122040190,122040251,122040273,14204026,14204057,142040108,142040115,142040239,142040285,142040289,142040293,162040156,162040178,162040208,162040210,162040217,162040323,162040355,162040372,162040395,17004043,162040421,17004096,170040109,14020154,50020150,6205078,70050124,82050125,90050310,94010242,94010267,10201080,102060204,114070105,150010327,17000046,170060159,194060155,230000138};

noc=Length[idc]

Extract the pubcat--.dat data  corresponding to the constant source ID's.

AbsoluteTiming[cst={};
Do[id=idc[[i]];iCon=Select[st,#[[1]]==id&];If[Length[iCon]>0,cst=Append[cst,iCon[[1]]]],{i,1,noc}
];ncst=Length[cst]]

Last[cst]

Here are the GB6 flux densities.

tcst=Transpose[cst];tcst[[13]]

Extract the corresponding daily flux densities and errors.

Read in Short Term (ST) daily flux values (pubfit10Jun16.dat)

dt=Import[ "e:/Lexar/capak/pubfit10Jun16.dat"];Length[dt]

Format of  "pubfit10Jun16.dat" data. 
1) Unique source serial number
2)  Scan number                     
3)  Daily flux density  (mJy)
4)  Daily flux density uncertainty  (mJy)        
5)  Modified Julian Day number (accurate to 4 decimal places)   MJD = JD - 2,400,000.5 
6)  Flag   (0= fit but (ave or wtd ave) flux too low compared to GB6, 1= successful fit)

Length[dt[[1]]]

dt[[1]]

dti=Select[dt,#[[1]]==idc[[20]]&&#[[5]]<46800.0&]

ndti=Length[dti]

tdti=Transpose[dti]

tdti[[3]]

tdti[[4]]

Here is the pubcat.dat data for source 20 with the GB 6 flux density and error in columns 13 & 14, respectively.

cst[[20]]

Here is the GB6 flux density.

tcst[[13]][[20]]

Here is the GB6 flux density error.

tcst[[14]][[20]]

dSE=(tdti[[3]]-tcst[[13]][[20]])

(Sqrt[((tdti[[4]])^2)+(((tcst[[14]][[1]]))^2)])

(tdti[[3]]-tcst[[13]][[20]])/(Sqrt[((tdti[[4]])^2)+((tcst[[14]][[20]])^2)])

lisDF86={};Do[dti=Select[dt,#[[1]]==idc[[i]]&&#[[5]]<46800.0&];ndti=Length[dti];If[ndti>0,tdti=Transpose[dti];dSE=(tdti[[3]]-tcst[[13]][[i]])/(Sqrt[((tdti[[4]])^2)+((tcst[[14]][[i]])^2)]);lisDF86=Append[lisDF86,dSE]],{i,1,noc}];Length[lisDF86]

nlisDF86=Length[Flatten[lisDF86]]

ph86=Histogram[Flatten[lisDF86],{-5,5,0.25},PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,ChartStyle->Gray,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(daily\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

pg86=Plot[nlisDF86*0.25*PDF[NormalDistribution[0,1],x],{x,-5,5},PlotStyle->{Black,Thick},Frame->True]

plHist86=Show[ph86,pg86,PlotRange->{{-5,5},All},Frame->True,FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(daily\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",12}],
		Style["Number",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Histogram ((((Subscript[S, daily]) - (Subscript[S, GB6])))/(Subscript[S, err])) for constant sources (epoch 1987)

Read list of constant source ID's (unresolved steep spectral index constant flux density sources).

idc={2040339,1404039,1404042,1404070,14040121,2204044,3004038,3004040,30040108,30040146,30040202,7404075,74040184,74040225,74040258,74040275,74040283,74040303,74040325,74040407,74040412,9404048,9404053,9404054,9404061,94040121,94040233,94040266,94040268,94040295,94040305,94040319,94040336,11404043,11404044,11404047,11404097,114040105,114040127,114040143,114040145,114040156,114040171,114040172,114040184,12204056,12204065,12204073,12204096,122040109,122040190,122040251,122040273,14204026,14204057,142040108,142040115,142040239,142040285,142040289,142040293,162040156,162040178,162040208,162040210,162040217,162040323,162040355,162040372,162040395,17004043,162040421,17004096,170040109,14020154,50020150,6205078,70050124,82050125,90050310,94010242,94010267,10201080,102060204,114070105,150010327,17000046,170060159,194060155,230000138};

noc=Length[idc]

Extract the pubcat--.dat data  corresponding to the constant source ID's.

AbsoluteTiming[cst={};
Do[id=idc[[i]];iCon=Select[st,#[[1]]==id&];If[Length[iCon]>0,cst=Append[cst,iCon[[1]]]],{i,1,noc}
];ncst=Length[cst]]

Last[cst]

Here are the GB6 flux densities.

tcst=Transpose[cst];tcst[[13]]

Extract the corresponding daily flux densities and errors.

Read in Short Term (ST) daily flux values (pubfit10Jun16.dat)

dt=Import[ "e:/Lexar/capak/pubfit10Jun16.dat"];Length[dt]

Format of  "pubfit10Jun16.dat" data. 
1) Unique source serial number
2)  Scan number                     
3)  Daily flux density  (mJy)
4)  Daily flux density uncertainty  (mJy)        
5)  Modified Julian Day number (accurate to 4 decimal places)   MJD = JD - 2,400,000.5 
6)  Flag   (0= fit but (ave or wtd ave) flux too low compared to GB6, 1= successful fit)

Length[dt[[1]]]

dt[[1]]

dti=Select[dt,#[[1]]==idc[[20]]&&#[[5]]>46800.0&]

ndti=Length[dti]

tdti=Transpose[dti]

tdti[[3]]

tdti[[4]]

Here is the pubcat.dat data for source 20 with the GB 6 flux density and error in columns 13 & 14, respectively.

cst[[20]]

Here is the GB6 flux density.

tcst[[13]][[20]]

Here is the GB6 flux density error.

tcst[[14]][[20]]

dSE=(tdti[[3]]-tcst[[13]][[20]])

(Sqrt[((tdti[[4]])^2)+(((tcst[[14]][[1]]))^2)])

(tdti[[3]]-tcst[[13]][[20]])/(Sqrt[((tdti[[4]])^2)+((tcst[[14]][[20]])^2)])

lisDF87={};Do[dti=Select[dt,#[[1]]==idc[[i]]&&#[[5]]>46800.0&];ndti=Length[dti];If[ndti>0,tdti=Transpose[dti];dSE=(tdti[[3]]-tcst[[13]][[i]])/(Sqrt[((tdti[[4]])^2)+((tcst[[14]][[i]])^2)]);lisDF87=Append[lisDF87,dSE]],{i,1,noc}];Length[lisDF87]

nlisDF87=Length[Flatten[lisDF87]]

ph87=Histogram[Flatten[lisDF87],{-5,5,0.25},PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,ChartStyle->Gray,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(daily\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

pg87=Plot[nlisDF87*0.25*PDF[NormalDistribution[0,1],x],{x,-5,5},PlotStyle->{Black,Thick},Frame->True]

plHist87=Show[ph87,pg87,PlotRange->{{-5,5},All},Frame->True,FrameLabel->{Style["\!\(\*FractionBox[\((\*SubscriptBox[\(S\), \(daily\)]\\  - \\ \*SubscriptBox[\(S\), \(GB6\)])\), SubscriptBox[\(S\), \(err\)]]\)",{"Helvetica-Bold",12}],
		Style["Number",{"Helvetica-Bold",12}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]]

Read in Short Term (ST) daily flux values (pubfit10Jun16.dat)  and examine S/N
(you need to edit the path for the location of the catalog on your computer)

dt=Import[ "e:/Lexar/capak/pubfit10Jun16.dat"];Length[dt]

Format of  "pubfit10Jun16.dat" data. 
1) Unique source serial number
2)  Scan number                     
3)  Daily flux density  (mJy)
4)  Daily flux density uncertainty  (mJy)        
5)  Modified Julian Day number (accurate to 4 decimal places)   MJD = JD - 2,400,000.5 
6)  Flag   (0= fit but (ave or wtd ave) flux too low compared to GB6, 1= successful fit)

Length[dt[[1]]]

tdt=Transpose[dt];r=tdt[[3]]/tdt[[4]];Median[r]

Mean[r]

ph=Histogram[r,PlotRange->All,Frame->True,Axes->False,AspectRatio->0.5,ChartStyle->Gray,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["((Subscript[S, daily])/(Subscript[S, err]))",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

(*Export["e:/Lexar/capak/examineVarLT_ST_19Jun16M7_HistogramDailyFluxToError.eps",ph,ImageSize->72*7];*)

dt[[1]]

Retrieve all ST source data with -(Subscript[Log, 10])(p-value) ≥ 3 (based on both epochs) and output source list with a header.

st[[1]]

ns

Column containg the 1986 epoch -(Subscript[Log, 10])(Subscript[(p-value), 86])

nP86=43;

Column containg the 1986 epoch -(Subscript[Log, 10])(Subscript[(p-value), 87])

nP87=44;

Set lower bound on - (Subscript[Log, 10])(p-value).

qMin=3;

Below msj is a list of the source parameters for the selected data.

msj={};nV=0;Do[If[st[[i,nP86]]≥qMin||st[[i,nP87]]≥qMin,nV++;msj=Append[msj,st[[i]]]],{i,1,ns}];nV

msj[[1]]

Last[msj]

Add an additional column to msj containing the larger of -(Subscript[Log, 10])(Subscript[(p-value), 86]) or -(Subscript[Log, 10])(Subscript[(p-value), 87]).

Do[q86=msj[[i,43]];q87=msj[[i,44]];qmax=If[q86>q87,q86,q87];msj[[i]]=Append[msj[[i]],qmax],{i,1,nV}];Last[msj]

Provide selected source list according to largest - Log(p-value).

msj=RotateLeft[Sort[RotateRight[msj,{0,1}],OrderedQ[{#2, #1}]&],{0,1}];

First[msj]

Choose particular source parameters to list in final table of selected sources.

printout=Table[Join[Take[msj[[i]],{1}],Take[msj[[i]],{9,10}],Take[msj[[i]],{21}],Take[msj[[i]],{22}],Take[msj[[i]],{13,14}],Take[msj[[i]],{23}],Take[msj[[i]],{30,31}],Take[msj[[i]],{24}],Take[msj[[i]],{32,33}],Take[msj[[i]],{40,41}],Take[msj[[i]],{42,45}]],{i,1,nV}];

First[printout]

Specify the header information for columns in the table of selected sources.

l1={{"ID","RA2000","DEC2000","B","L","SGB6","Serr","no86","Swt86","Swt86err","no87","Swt87","Swt87err","MI86","MI87","-(Subscript[Log, 10])(p-value)","-(Subscript[Log, 10])(p-value)_86","-(Subscript[Log, 10])(p-value)_87","-(Subscript[Log, 10])(p-value)_max"},{" "," "," "," "," ","mJy","mJy","","mJy","mJy"," ","mJy","mJy"," corrected","corrected "," "," "," "," "}}

Subset of top 50 ST variables selected sources.

smj=MatrixForm[Join[l1,Take[Transpose[Transpose[printout]],50]]]

Create the full table of the selected sources properties with header information.

data=Join[l1,printout];

To export the data file remove the (*  *) at the ends of the next command.

(*Export["e:/Lexar/capak/top_LogPvalue_selected.dat",data,"Table"];*)

Plot daily flux densities for sources with highest -(Subscript[Log, 10])(p-value)  (based on both epochs) and save plots.

sorc=Drop[Import["e:/Lexar/capak/top_LogPvalue_selected.dat"],3];

FullForm[sorc[[1]]]

l1={{"ID","RA2000","DEC2000","B","L","SGB6","Serr","no86","Swt86","Swt86err","no87","Swt87","Swt87err","MI86","MI87","-(Subscript[Log, 10])(p-value)","-(Subscript[Log, 10])(p-value)_86","-(Subscript[Log, 10])(p-value)_87","-(Subscript[Log, 10])(p-value)_max"},{" "," "," "," "," ","mJy","mJy"," ","mJy","mJy"," ","mJy","mJy"," corrected","corrected "," "," "," "," "}}

nos=Length[sorc]

Read in Short Term (ST) daily flux values (pubfit10Jun16.dat)

data=Import[ "e:/Lexar/capak/pubfit10Jun16.dat"];
nd=Length[data]

Format of  "pubfit10Jun16.dat" data. 
1) Unique source serial number
2)  Scan number                     
3)  Daily flux density  (mJy)
4)  Daily flux density uncertainty  (mJy)        
5)  Modified Julian Day number (accurate to 4 decimal places)   MJD = JD - 2,400,000.5 
6)  Flag   (0= fit but (ave or wtd ave) flux too low compared to GB6, 1= successful fit)

Length[data[[1]]]

FullForm[data[[1]]]

FullForm[sorc[[1]][[1]]]

sorc[[5]][[1]]

nsource=4;k=1;j=0;noplot=nsource*2;pl=Table[1,{i,noplot}];mjd1=2446740.5-2400000.5;mjd2=2447067.5-2400000.5;
While[k≤nsource,n1=0;n2=0;
name=sorc[[k]][[1]];ra=sorc[[k]][[2]];dec=sorc[[k]][[3]];fl=Select[data,#[[1]]==name&];nfl=Length[fl];
If[nfl≥2,tfl=Transpose[fl];upper=Max[tfl[[3]]];flA=Sort[Transpose[{tfl[[5]],tfl[[3]],tfl[[4]]}],#1[[1]]<#2[[1]]&];Do[If[flA[[i]][[1]]≥46800,n2++,n1++],{i,1,nfl}
]
];
If[n1>0,j++;fl1=Table[{{flA[[i]][[1]]-mjd1,flA[[i]][[2]]},ErrorBar[flA[[i]][[3]]]},{i,1,n1}];pl[[j]]=
ErrorListPlot[fl1,Frame->True,PlotRange->{{-1,36},{-0.02,1.15*upper+0.02}},Axes->False,AspectRatio->0.5,PlotStyle->Black,PlotLabel->Style[ToString[name]<>"   (Subscript[RA, 2000]) = "<>ra<>"   DEC = "<>dec<>"   k = "<>ToString[k]<>"\n NLpvalue86="<>ToString[sorc[[k]][[17]]]<>"  NLpvalue87="<>ToString[sorc[[k]][[18]]]<>"  b="<>ToString[sorc[[k]][[4]]]<>"  l="<>ToString[sorc[[k]][[5]]],{"Helvetica-Bold",8}],FrameLabel->{Style["Time from start of epoch (MJD-"<>ToString[mjd1]<>")",{"Helvetica-Bold",13}],Style["Flux density  (mJy)",{"Helvetica-Bold",13}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]
]
];If[n2>0,j++;fl2=Table[{{flA[[i]][[1]]-mjd2,flA[[i]][[2]]},ErrorBar[flA[[i]][[3]]]},{i,n1+1,n1+n2}];
pl[[j]]=
ErrorListPlot[fl2,Frame->True,PlotRange->{{-1,36},{-0.02,1.15*upper+0.02}},Axes->False,AspectRatio->0.5,PlotStyle->Black,PlotLabel->Style[ToString[name]<>"   (Subscript[RA, 2000]) = "<>ra<>"   DEC = "<>dec<>"   k = "<>ToString[k]<>"\n NLpvalue86="<>ToString[sorc[[k]][[17]]]<>"  NLpvalue87="<>ToString[sorc[[k]][[18]]]<>"  b="<>ToString[sorc[[k]][[4]]]<>"  l="<>ToString[sorc[[k]][[5]]],{"Helvetica-Bold",8}],FrameLabel->{Style["Time from start of epoch (MJD-"<>ToString[mjd2]<>")",{"Helvetica-Bold",13}],Style["Flux density  (mJy)",{"Helvetica-Bold",13}]},BaseStyle->{"Helvetica-Bold",13},FrameStyle->Thickness[0.006]
]
];k++
]

If[j>0,plAll=Show[GraphicsGrid[If[EvenQ[j],Table[{pl[[i]],pl[[i+1]]},{i,1,j,2}],Join[Table[{pl[[i]],pl[[i+1]]},{i,1,IntegerPart[j/2]*2,2}],{{pl[[IntegerPart[j/2]*2+1]]}}]],Spacings->{Scaled[-0.25],Scaled[0.1]}]]]

(*Export["e:/Lexar/capak/top_LogPvalue_selected.eps",plAll,ImageSize->72*11];*)

Plot a Parabolic Projection In Celestial Coordinates  (1950) -(Subscript[Log, 10])(p-value) < 0.1 (1986 epoch)

RA defined positive measured eastward from the vernal equinox

k=(π/180);

xc[ϕc_,θc_]:=-ϕc (2 Cos[((2 k θc)/3)] - 1)

yc[ϕc_,θc_]:=180 Sin[((k θc)/3)]

N[xc[180,0]]

l1=ParametricPlot[{xc[-180,θc],yc[-180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l2=ParametricPlot[{xc[-120,θc],yc[-120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l3=ParametricPlot[{xc[-60,θc],yc[-60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l4=ParametricPlot[{xc[60,θc],yc[60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l5=ParametricPlot[{xc[120,θc],yc[120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l6=ParametricPlot[{xc[180,θc],yc[180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];

b1=ParametricPlot[{xc[ϕc,60],yc[ϕc,60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b2=ParametricPlot[{xc[ϕc,30],yc[ϕc,30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b3=ParametricPlot[{xc[ϕc,-30],yc[ϕc,-30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b4=ParametricPlot[{xc[ϕc,-60],yc[ϕc,-60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];

Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4}]

ns

Format of "pubcat10Jun16.dat" that contains summary information and statistics on the daily flux density measurements.
1)   Serial Number    014050196
2)    Δθ = (daily scan data derived position - GB6 position)  arcmin                     
3)   HH:MM:SS.S        GB6 Right Ascension (RA) (B1950)
4)   SS.S                     Uncertainty in right ascension (s)
5)   SDD:MM:SS         GB6 Declination (Dec) (B1950)
6)   SS                        Uncertainty in declination (arcsec)
7)    GB6 RA      (B1950, decimal degrees)
8)    GB6 Dec    (B1950, decimal degrees)
9)     HH:MM:SS.S        GB6 RA   (J2000) 
10)   SDD:MM:SS.SS  GB6 Dec  (J2000) 
11)  GB6 RA     (J2000, decimal degrees)
12)  GB6 Dec   (J2000, decimal degrees)        
13)  GB6 Flux density            (mJy)
14)  GB6 Flux density Error    (mJy)
15)  Number of confusing sources subtracted
16)  Maximum Scan Flux density    (mJy)
17)  Fitted RA in minutes of time        (1950)
18)  Fitted RA in minutes of time error 
19)  Fitted Dec in decimal degrees    (1950)
20)  Fitted Dec in decimal degrees error  
21)  Galactic Latitude       (decimal degrees)
22)  Galactic Longitude    (decimal degrees)
23)  Number of fits in epoc 1986
24)  Number of fits in epoc 1987
25)  Average Flux density (both epochs)      (mJy)
26)  Average Flux density epoc 1986            (mJy)
27)  Average Flux density epoc 1987            (mJy)
28)  Weighted average Flux density (both epochs)    (mJy)
29)    Error in weighted average                                  (mJy)
30)  Weighted average Flux density epoc 1986          (mJy)
31)    Error in weighted average1986                          (mJy)
32)  Weighted average Flux density epoc 1987          (mJy)
33)    Error in weighted average 1987                         (mJy)
34)  Modulation Index  epoc 1986  given by (1/(OverBar[(Subscript[S, wt86])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i86])-(OverBar[(Subscript[S, wt86])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))
35)  Modulation Index  epoc 1987 given by (1/(OverBar[(Subscript[S, wt87])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i87])-(OverBar[(Subscript[S, wt87])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))  
36)  Modulation Index due to noise from 1000 simulations for epoc 1986   (-9 means ignore)
37)  Modulation Index due to noise from 1000 simulations for epoc 1987   (-9 means ignore)
38)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1986   (-9 means ignore)
39)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1987   (-9 means ignore) 
40)  Corrected Modulation Index epoc 1986
41)  Corrected Modulation Index epoc 1987
42)  -(Subscript[log, 10]) (p-value)  
43)  -(Subscript[log, 10]) (p-value)  for epoc  1   
44)  -(Subscript[log, 10]) (p-value) for epoc 2

First entry in the short term variability data set (pubcat10Jun16.dat)

st[[1]]

Histogram of the RA values in HH.HHHH of the selected data.

plRAHist=Histogram[Transpose[st][[11]]/15,{0.5}]

Histogram of the RA values in DD.DDDD of the selected data.

plRAHist=Histogram[Transpose[st][[11]],{5.0},PlotRange->{{0,360},All},Frame->True,ChartStyle->EdgeForm[],Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["RA deg (J2000)",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Set lower bound on desired -(Subscript[Log, 10])(p-value).

qMin=0.1;

Specify column for desired -(Subscript[Log, 10])(p-value).

colLogQ=43;

Next we select the sources with desired -(Subscript[Log, 10])(p-value).

nsel=0;vg={};Do[If[st[[i]][[colLogQ]]<qMin,nsel++;vg=Append[vg,{st[[i]][[7]],st[[i]][[8]]}],],{i,1,ns}];nsel

Construct a table of RA, DEC pairs for each selected source

y=Table[{xc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]],yc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]]},{i,1,nsel}];

Plot the locations of the selected sources in RA, DEC

p1=ListPlot[y,PlotRange->{{-180,180},{-90,90}}];

radecpl=Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4,p1},DisplayFunction->$DisplayFunction]

Plot a Parabolic Projection In Celestial Coordinates  (1950) -(Subscript[Log, 10])(p-value) ≥2 (1987 epoch).

RA defined positive measured eastward from the vernal equinox

k=(π/180);

xc[ϕc_,θc_]:=-ϕc (2 Cos[((2 k θc)/3)] - 1)

yc[ϕc_,θc_]:=180 Sin[((k θc)/3)]

N[xc[180,0]]

l1=ParametricPlot[{xc[-180,θc],yc[-180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l2=ParametricPlot[{xc[-120,θc],yc[-120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l3=ParametricPlot[{xc[-60,θc],yc[-60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l4=ParametricPlot[{xc[60,θc],yc[60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l5=ParametricPlot[{xc[120,θc],yc[120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l6=ParametricPlot[{xc[180,θc],yc[180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];

b1=ParametricPlot[{xc[ϕc,60],yc[ϕc,60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b2=ParametricPlot[{xc[ϕc,30],yc[ϕc,30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b3=ParametricPlot[{xc[ϕc,-30],yc[ϕc,-30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b4=ParametricPlot[{xc[ϕc,-60],yc[ϕc,-60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];

Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4}]

ns

First entry in the short term variability data set (pubcat10Jun16.dat)

st[[1]]

Format of "pubcat10Jun16.dat" that contains summary information and statistics on the daily flux density measurements.
1)   Serial Number    014050196
2)    Δθ = (daily scan data derived position - GB6 position)  arcmin                     
3)   HH:MM:SS.S        GB6 Right Ascension (RA) (B1950)
4)   SS.S                     Uncertainty in right ascension (s)
5)   SDD:MM:SS         GB6 Declination (Dec) (B1950)
6)   SS                        Uncertainty in declination (arcsec)
7)    GB6 RA      (B1950, decimal degrees)
8)    GB6 Dec    (B1950, decimal degrees)
9)     HH:MM:SS.S        GB6 RA   (J2000) 
10)   SDD:MM:SS.SS  GB6 Dec  (J2000) 
11)  GB6 RA     (J2000, decimal degrees)
12)  GB6 Dec   (J2000, decimal degrees)        
13)  GB6 Flux density            (mJy)
14)  GB6 Flux density Error    (mJy)
15)  Number of confusing sources subtracted
16)  Maximum Scan Flux density    (mJy)
17)  Fitted RA in minutes of time        (1950)
18)  Fitted RA in minutes of time error 
19)  Fitted Dec in decimal degrees    (1950)
20)  Fitted Dec in decimal degrees error  
21)  Galactic Latitude       (decimal degrees)
22)  Galactic Longitude    (decimal degrees)
23)  Number of fits in epoc 1986
24)  Number of fits in epoc 1987
25)  Average Flux density (both epochs)      (mJy)
26)  Average Flux density epoc 1986            (mJy)
27)  Average Flux density epoc 1987            (mJy)
28)  Weighted average Flux density (both epochs)    (mJy)
29)    Error in weighted average                                  (mJy)
30)  Weighted average Flux density epoc 1986          (mJy)
31)    Error in weighted average1986                          (mJy)
32)  Weighted average Flux density epoc 1987          (mJy)
33)    Error in weighted average 1987                         (mJy)
34)  Modulation Index  epoc 1986  given by (1/(OverBar[(Subscript[S, wt86])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i86])-(OverBar[(Subscript[S, wt86])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))
35)  Modulation Index  epoc 1987 given by (1/(OverBar[(Subscript[S, wt87])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i87])-(OverBar[(Subscript[S, wt87])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))  
36)  Modulation Index due to noise from 1000 simulations for epoc 1986   (-9 means ignore)
37)  Modulation Index due to noise from 1000 simulations for epoc 1987   (-9 means ignore)
38)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1986   (-9 means ignore)
39)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1987   (-9 means ignore) 
40)  Corrected Modulation Index epoc 1986
41)  Corrected Modulation Index epoc 1987
42)  -(Subscript[log, 10]) (p-value)  
43)  -(Subscript[log, 10]) (p-value)  for epoc  1   
44)  -(Subscript[log, 10]) (p-value) for epoc 2

Set lower bound on desired -(Subscript[Log, 10])(p-value).

qMin=2;

Specify column for desired -(Subscript[Log, 10])(p-value).

colLogQ=44;

Next we select the sources with desired -(Subscript[Log, 10])(p-value).

nsel=0;vg={};Do[If[st[[i]][[colLogQ]]≥qMin,nsel++;vg=Append[vg,{st[[i]][[7]],st[[i]][[8]]}],],{i,1,ns}];nsel

Construct a table of RA, DEC pairs for each selected source

y=Table[{xc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]],yc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]]},{i,1,nsel}];

Plot the locations of the selected sources in RA, DEC

p1=ListPlot[y,PlotRange->{{-180,180},{-90,90}}];

radecpl=Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4,p1},DisplayFunction->$DisplayFunction]

Plot a Parabolic Projection In Galactic Coordinates  (1950) -(Subscript[Log, 10])(p-value) ≥1 (based on 1986 or 1987 epochs).

RA defined positive measured eastward from the vernal equinox

k=(π/180);

xc[ϕc_,θc_]:=-ϕc (2 Cos[((2 k θc)/3)] - 1)

yc[ϕc_,θc_]:=180 Sin[((k θc)/3)]

N[xc[180,0]]

l1=ParametricPlot[{xc[-180,θc],yc[-180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l2=ParametricPlot[{xc[-120,θc],yc[-120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l3=ParametricPlot[{xc[-60,θc],yc[-60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l4=ParametricPlot[{xc[60,θc],yc[60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l5=ParametricPlot[{xc[120,θc],yc[120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l6=ParametricPlot[{xc[180,θc],yc[180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];

b1=ParametricPlot[{xc[ϕc,60],yc[ϕc,60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b2=ParametricPlot[{xc[ϕc,30],yc[ϕc,30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b3=ParametricPlot[{xc[ϕc,-30],yc[ϕc,-30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b4=ParametricPlot[{xc[ϕc,-60],yc[ϕc,-60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];

Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4}]

ns

First entry in the short term variability data set (pubcat10Jun16.dat)

st[[1]]

Format of "pubcat10Jun16.dat" that contains summary information and statistics on the daily flux density measurements.
1)   Serial Number    014050196
2)    Δθ = (daily scan data derived position - GB6 position)  arcmin                     
3)   HH:MM:SS.S        GB6 Right Ascension (RA) (B1950)
4)   SS.S                     Uncertainty in right ascension (s)
5)   SDD:MM:SS         GB6 Declination (Dec) (B1950)
6)   SS                        Uncertainty in declination (arcsec)
7)    GB6 RA      (B1950, decimal degrees)
8)    GB6 Dec    (B1950, decimal degrees)
9)     HH:MM:SS.S        GB6 RA   (J2000) 
10)   SDD:MM:SS.SS  GB6 Dec  (J2000) 
11)  GB6 RA     (J2000, decimal degrees)
12)  GB6 Dec   (J2000, decimal degrees)        
13)  GB6 Flux density            (mJy)
14)  GB6 Flux density Error    (mJy)
15)  Number of confusing sources subtracted
16)  Maximum Scan Flux density    (mJy)
17)  Fitted RA in minutes of time        (1950)
18)  Fitted RA in minutes of time error 
19)  Fitted Dec in decimal degrees    (1950)
20)  Fitted Dec in decimal degrees error  
21)  Galactic Latitude       (decimal degrees)
22)  Galactic Longitude    (decimal degrees)
23)  Number of fits in epoc 1986
24)  Number of fits in epoc 1987
25)  Average Flux density (both epochs)      (mJy)
26)  Average Flux density epoc 1986            (mJy)
27)  Average Flux density epoc 1987            (mJy)
28)  Weighted average Flux density (both epochs)    (mJy)
29)    Error in weighted average                                  (mJy)
30)  Weighted average Flux density epoc 1986          (mJy)
31)    Error in weighted average1986                          (mJy)
32)  Weighted average Flux density epoc 1987          (mJy)
33)    Error in weighted average 1987                         (mJy)
34)  Modulation Index  epoc 1986  given by (1/(OverBar[(Subscript[S, wt86])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i86])-(OverBar[(Subscript[S, wt86])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))
35)  Modulation Index  epoc 1987 given by (1/(OverBar[(Subscript[S, wt87])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i87])-(OverBar[(Subscript[S, wt87])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))  
36)  Modulation Index due to noise from 1000 simulations for epoc 1986   (-9 means ignore)
37)  Modulation Index due to noise from 1000 simulations for epoc 1987   (-9 means ignore)
38)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1986   (-9 means ignore)
39)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1987   (-9 means ignore) 
40)  Corrected Modulation Index epoc 1986
41)  Corrected Modulation Index epoc 1987
42)  -(Subscript[log, 10]) (p-value)  
43)  -(Subscript[log, 10]) (p-value)  for epoc  1   
44)  -(Subscript[log, 10]) (p-value) for epoc 2

Set lower bound on desired -(Subscript[Log, 10])(p-value).

qMin=1;

Specify column for desired -(Subscript[Log, 10])(p-value).

colLogQ86=43;colLogQ87=44;

Next we select the sources with desired -(Subscript[Log, 10])(p-value).

nsel=0;vg={};Do[If[st[[i]][[colLogQ86]]≥qMin||st[[i]][[colLogQ87]]≥qMin,nsel++;vg=Append[vg,{st[[i]][[22]],st[[i]][[21]]}]],{i,1,ns}];nsel

Construct a table of galactic coordinates (l, b) pairs for each selected source

y=Table[{xc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]],yc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]]},{i,1,nsel}];

Plot the locations of the selected sources in l, b

p1=ListPlot[y,PlotRange->{{-180,180},{-90,90}}];

Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4,p1}]

Plot a Parabolic Projection In Celestial Coordinates  (1950) for (n86+n87) ≥ 12.

RA defined positive measured eastward from the vernal equinox

k=(π/180);

xc[ϕc_,θc_]:=-ϕc (2 Cos[((2 k θc)/3)] - 1)

yc[ϕc_,θc_]:=180 Sin[((k θc)/3)]

N[xc[180,0]]

l1=ParametricPlot[{xc[-180,θc],yc[-180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l2=ParametricPlot[{xc[-120,θc],yc[-120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l3=ParametricPlot[{xc[-60,θc],yc[-60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l4=ParametricPlot[{xc[60,θc],yc[60,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l5=ParametricPlot[{xc[120,θc],yc[120,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];l6=ParametricPlot[{xc[180,θc],yc[180,θc]},{θc,-90,90},PlotRange->{{-180,180},{-90,90}}];

b1=ParametricPlot[{xc[ϕc,60],yc[ϕc,60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b2=ParametricPlot[{xc[ϕc,30],yc[ϕc,30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b3=ParametricPlot[{xc[ϕc,-30],yc[ϕc,-30]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];b4=ParametricPlot[{xc[ϕc,-60],yc[ϕc,-60]},{ϕc,-180,180},PlotRange->{{-180,180},{-90,90}}];

Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4}]

ns

First entry in the short term variability data set (pubcat10Jun16.dat)

st[[1]]

Format of "pubcat10Jun16.dat" that contains summary information and statistics on the daily flux density measurements.
1)   Serial Number    014050196
2)    Δθ = (daily scan data derived position - GB6 position)  arcmin                     
3)   HH:MM:SS.S        GB6 Right Ascension (RA) (B1950)
4)   SS.S                     Uncertainty in right ascension (s)
5)   SDD:MM:SS         GB6 Declination (Dec) (B1950)
6)   SS                        Uncertainty in declination (arcsec)
7)    GB6 RA      (B1950, decimal degrees)
8)    GB6 Dec    (B1950, decimal degrees)
9)     HH:MM:SS.S        GB6 RA   (J2000) 
10)   SDD:MM:SS.SS  GB6 Dec  (J2000) 
11)  GB6 RA     (J2000, decimal degrees)
12)  GB6 Dec   (J2000, decimal degrees)        
13)  GB6 Flux density            (mJy)
14)  GB6 Flux density Error    (mJy)
15)  Number of confusing sources subtracted
16)  Maximum Scan Flux density    (mJy)
17)  Fitted RA in minutes of time        (1950)
18)  Fitted RA in minutes of time error 
19)  Fitted Dec in decimal degrees    (1950)
20)  Fitted Dec in decimal degrees error  
21)  Galactic Latitude       (decimal degrees)
22)  Galactic Longitude    (decimal degrees)
23)  Number of fits in epoc 1986
24)  Number of fits in epoc 1987
25)  Average Flux density (both epochs)      (mJy)
26)  Average Flux density epoc 1986            (mJy)
27)  Average Flux density epoc 1987            (mJy)
28)  Weighted average Flux density (both epochs)    (mJy)
29)    Error in weighted average                                  (mJy)
30)  Weighted average Flux density epoc 1986          (mJy)
31)    Error in weighted average1986                          (mJy)
32)  Weighted average Flux density epoc 1987          (mJy)
33)    Error in weighted average 1987                         (mJy)
34)  Modulation Index  epoc 1986  given by (1/(OverBar[(Subscript[S, wt86])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i86])-(OverBar[(Subscript[S, wt86])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))
35)  Modulation Index  epoc 1987 given by (1/(OverBar[(Subscript[S, wt87])]))(Sqrt[((∑(((Subscript[w, i])((Subscript[S, i87])-(OverBar[(Subscript[S, wt87])])))^2))/(∑(Subscript[w, i])))]), where (Subscript[w, i]) = (1/((Subscript[σ, i])^2))  
36)  Modulation Index due to noise from 1000 simulations for epoc 1986   (-9 means ignore)
37)  Modulation Index due to noise from 1000 simulations for epoc 1987   (-9 means ignore)
38)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1986   (-9 means ignore)
39)  Std deviation from 1000 simulations for Modulation index due to noise for epoc 1987   (-9 means ignore) 
40)  Corrected Modulation Index epoc 1986
41)  Corrected Modulation Index epoc 1987
42)  -(Subscript[log, 10]) (p-value)  
43)  -(Subscript[log, 10]) (p-value)  for epoc  1   
44)  -(Subscript[log, 10]) (p-value) for epoc 2

Next we select the desired sources.

nsel=0;vg={};Do[If[(st[[i]][[23]]+st[[i]][[23]])≥12.0,nsel++;vg=Append[vg,{st[[i]][[7]],st[[i]][[8]]}],],{i,1,ns}];nsel

Construct a table of RA, DEC pairs for each selected source

y=Table[{xc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]],yc[If[vg[[i]][[1]]≤180,vg[[i]][[1]],(vg[[i]][[1]]-360)],vg[[i]][[2]]]},{i,1,nsel}];

Plot the locations of the selected sources in RA, DEC

p1=ListPlot[y,PlotRange->{{-180,180},{-90,90}}];

radecpl=Show[{l1,l2,l3,l4,l5,l6,b1,b2,b3,b4,p1},DisplayFunction->$DisplayFunction]