Commits
Takeshi Nakazato authored and Ville Suoranta committed 2af391c38ea Merge
32 32 | ScalarFormatter) |
33 33 | from six.moves import input, range |
34 34 | |
35 35 | # CAS-13722, CAS-13385 |
36 36 | import warnings |
37 37 | import matplotlib.cbook |
38 38 | #warnings.filterwarnings("ignore",category=matplotlib.cbook.MatplotlibDeprecationWarning) |
39 39 | |
40 40 | ##------------------------------------------------------------------------------------------------------ |
41 41 | ##--- Increment the micro version number with each update to task_plotbandpass.py, adjust the major --- |
42 - | ##--- and minor version numbers to match the analusisUtils version that is last synced with... --- |
42 + | ##--- and minor version numbers to match the analysisUtils version that is last synced with... --- |
43 43 | ##------------------------------------------------------------------------------------------------------ |
44 - | TASK_PLOTBANDPASS_REVISION_STRING = "2.3.0" |
44 + | TASK_PLOTBANDPASS_REVISION_STRING = "2.20.0" #### incremented from 2.3.0 in 2024Aug |
45 45 | ##------------------------------------------------------------------------------------------------------ |
46 46 | ##--- Update this version string whenever task_plotbandpass.py is synced with plotbandpass3.py from --- |
47 47 | ##--- Todd's analysisUtils. This will allow for tracking efforts to keep the versions synced. --- |
48 48 | ##------------------------------------------------------------------------------------------------------ |
49 - | PLOTBANDPASS_REVISION_STRING = "$Id: plotbandpass3.py,v 2.3 2023/01/06 17:41:16 thunter Exp $" |
49 + | PLOTBANDPASS_REVISION_STRING = "Id: plotbandpass3.py,v 2.20 2024/09/04 16:00:03 thunter Exp" |
50 + | |
50 51 | TOP_MARGIN = 0.25 # Used if showatm=T or showtksy=T |
51 52 | BOTTOM_MARGIN = 0.25 # Used if showfdm=T |
52 53 | MAX_ATM_CALC_CHANNELS = 512 |
53 54 | |
54 55 | markeredgewidth = 0.0 |
56 + | maxAltitude = 60 # for ozone, in km, this is the default value of the parameter in the analysisUtils version |
55 57 | |
56 58 | # This is a color sequence found online which has distinguishable colors |
57 59 | overlayColorsSequence = [ |
58 60 | [0.00, 0.00, 0.00], |
59 61 | [0.00, 0.00, 1.00], |
60 62 | [0.00, 0.50, 0.00], |
61 63 | [1.00, 0.00, 0.00], |
62 64 | [0.00, 0.75, 0.75], |
63 65 | # [0.75, 0.00, 0.75], # magenta, same as atmcolor |
64 66 | [0.25, 0.25, 0.25], |
69 71 | [0.00, 1.00, 0.00], |
70 72 | [0.76, 0.57, 0.17], |
71 73 | [0.54, 0.63, 0.22], |
72 74 | [0.34, 0.57, 0.92], |
73 75 | [1.00, 0.10, 0.60], |
74 76 | # [0.88, 0.75, 0.73], invisible |
75 77 | [0.10, 0.49, 0.47], |
76 78 | [0.66, 0.34, 0.65], |
77 79 | [0.99, 0.41, 0.23]] |
78 80 | overlayColorsList = overlayColorsSequence.copy() |
79 - | overlayColorsList += overlayColorsList + overlayColorsList # 17*3 = 51 total colors |
80 - | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time |
81 - | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time |
82 - | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time |
83 - | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time |
81 + | overlayColorsList += overlayColorsList + overlayColorsList # 17*3 = 34 total color entries |
82 + | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time 34*3 =102 |
83 + | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time 102*3=306 |
84 + | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time 306*3=918 |
85 + | overlayColorsList += overlayColorsList + overlayColorsList # try to support antenna,time 918*3=2754 entries |
84 86 | |
85 87 | # Enumeration to keep track of plot pages |
86 88 | PAGE_ANT = 0 |
87 89 | PAGE_SPW = 1 |
88 90 | PAGE_TIME = 2 |
89 91 | PAGE_AP = 3 |
90 92 | |
91 93 | # Used to parse command line arguments |
92 94 | myValidCharacterList = ['~', ',', ' ', '*',] + [str(m) for m in range(10)] |
93 95 | myValidCharacterListWithBang = ['~', ',', ' ', '*', '!',] + [str(m) for m in range(10)] |
94 96 | LARGE_POSITIVE = +1e20 |
95 97 | LARGE_NEGATIVE = -1e20 |
96 98 | maxAntennaNamesAcrossTheTop = 17 |
97 99 | maxTimesAcrossTheTop = 13 # 17 for HH:MM, reduced by 1 below for subplot=11 |
98 100 | antennaVerticalSpacing = 0.018 # 0.016 |
99 101 | antennaHorizontalSpacing = 0.05 |
100 102 | xstartTitle = 0.07 |
101 103 | ystartTitle = 0.955 |
102 104 | xstartPolLabel = 0.05 |
103 - | ystartOverlayLegend = 0.933 |
105 + | ystartOverlayLegend = 0.931 |
104 106 | opaqueSky = 270. # Kelvin, used for scaling TebbSky |
105 107 | |
106 108 | developerEmail = "thunter@nrao.edu" |
107 109 | |
108 110 | #class Polarization: |
109 111 | # taken from Stokes.h in casa, for reference only |
110 112 | # (Undefined, I,Q,U,V,RR,RL,LR,LL,XX,XY,YX,YY) = range(13) |
111 113 | |
112 114 | def version(showfile=True): |
113 115 | """ |
496 498 | mytb.close() |
497 499 | donetime = time.time() |
498 500 | return(nPolarizations) |
499 501 | |
500 502 | def getnspw(mymsmd): |
501 503 | # if (casadef.subversion_revision > '22653'): |
502 504 | return(mymsmd.nspw(False)) |
503 505 | # else: |
504 506 | # return(mymsmd.nspw()) |
505 507 | |
506 - | |
507 - | def drawOverlayTimeLegends(xframe,firstFrame,xstartTitle,ystartTitle,caltable,titlesize, |
508 - | fieldIndicesToPlot,ispwInCalTable,uniqueTimesPerFieldPerSpw, |
509 - | timerangeListTimes, solutionTimeThresholdSeconds,debugSloppyMatch, |
510 - | ystartOverlayLegend,debug,mysize, fieldsToPlot,myUniqueColor, |
508 + | def drawOverlayTimeLegends(xframe, firstFrame, xstartTitle, ystartTitle, caltable, titlesize, |
509 + | fieldIndicesToPlot, ispwInCalTable, uniqueTimesPerFieldPerSpw, |
510 + | timerangeListTimes, solutionTimeThresholdSeconds, debugSloppyMatch, |
511 + | ystartOverlayLegend, debug, mysize, fieldsToPlot, myUniqueColor, |
511 512 | timeHorizontalSpacing, fieldIndex, overlayColors, |
512 513 | antennaVerticalSpacing, overlayAntennas, |
513 - | timerangeList, caltableTitle, |
514 - | mytime, scansToPlot, scansForUniqueTimes): |
514 + | timerangeList, caltableTitle, mytime, scansToPlot, |
515 + | scansForUniqueTimes, uniqueSpwsInCalTable, uniqueTimes): |
515 516 | """ |
516 517 | Draws the legend at the top of the page, if it is the correct time to do so, |
517 518 | including the overlayTimes, the 'UT' label, and the caltable name. |
519 + | xframe: an integer as a subplot specifier, like 223 for the third panel of a 2x2 plot |
520 + | firstFrame: an integer as a subplot specifier, like 221 for the first panel of a 2x2 plot |
521 + | caltable: name of caltable on disk |
522 + | titlesize: font size as an integer string, e.g. '7' |
523 + | fieldIndicesToPlot: a list of integer indexes to the parent list fieldsToPlot, e.g. [0] |
524 + | fieldsToPlot: list of measurement set field IDs to plot, this can be different (shorter) than parent list fieldsToPlot |
525 + | scansToPlot: a list of measurement set scan numbers from which solutions should be plotted |
526 + | scansForUniqueTimes: a list of scans corresponding to the parent list of uniqueTimes |
527 + | ispwInCalTable: an integer index for the single desired spw in the spw list uniqueSpwsInCalTable, e.g. 0 |
528 + | uniqueTimesPerFieldPerSpw: a list of list of lists corresponding to [spwIndex][fieldIndex][floating-point times...] |
529 + | uniqueTimes: a list of all unique floating-point times in the caltable |
530 + | timerangeListTimes: a list of unique values of floating point MJD seconds that were requested to be plotted |
531 + | timerangeList: a list of timerange indices that were requested to be plotted |
532 + | uniqueSpwsInCalTable: a list of spw IDs, but only used for print statements |
533 + | mytime: the index of the timerange in the list of uniqueTimes in the caltable that was |
534 + | being examined when this function was called in the parent 'while loop' over uniqueTimes |
518 535 | """ |
519 536 | # debugSloppyMatch=True |
537 + | myspw = uniqueSpwsInCalTable[ispwInCalTable] |
538 + | if debug: |
539 + | print("len(timerangeList)=%d, len(timerangeListTimes)=%d" % (len(timerangeList), len(timerangeListTimes))) |
540 + | print("timerangeListTimes: %s" % (['%.3f'%(i) for i in timerangeListTimes])) |
541 + | print("xframe=%d, firstFrame=%d, ispwInCalTable=%d, spw=%d, fieldIndicesToPlot=%s, fieldsToPlot=%s" % (xframe,firstFrame,ispwInCalTable,myspw,fieldIndicesToPlot,fieldsToPlot)) |
542 + | print("uniqueTimesPerFieldPerSpw[%d][%d] = %s" % (ispwInCalTable,fieldIndicesToPlot[0],uniqueTimesPerFieldPerSpw[ispwInCalTable][fieldIndicesToPlot[0]])) |
520 543 | if (xframe == firstFrame): |
521 544 | # draw title including caltable name |
522 545 | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, color='k', transform=pb.gcf().transFigure) |
523 546 | # support multi-fields with overlay='time' |
524 547 | uTPFPS = [] |
525 548 | uTPFPStimerange = [] |
526 549 | |
527 - | # Find all timerange integers for all fields, not just the ones that were plotted |
550 + | # Find all timerange indices for all fields, not just the ones that were plotted |
528 551 | allTimeranges = [] |
529 552 | for f in range(len(uniqueTimesPerFieldPerSpw[ispwInCalTable])): |
530 553 | for t in uniqueTimesPerFieldPerSpw[ispwInCalTable][f]: |
531 - | if (t in timerangeListTimes): |
532 - | allTimeranges.append(list(timerangeListTimes).index(t)) |
533 - | |
554 + | if (np.min(np.abs(timerangeListTimes-t)) < solutionTimeThresholdSeconds): ### added 2024Aug |
555 + | closestIndex = np.argmin(np.abs(timerangeListTimes-t)) ### added 2024Aug |
556 + | allTimeranges.append(closestIndex) ### added 2024Aug |
557 + | |
534 558 | allTimeranges = list(np.sort(np.unique(allTimeranges))) |
535 - | |
559 + | # allTimeranges is a list of integers. |
560 + | # The length of allTimeranges will generally be shorter than timerangeListTimes. |
561 + | |
536 562 | for f in fieldIndicesToPlot: |
537 - | for t in uniqueTimesPerFieldPerSpw[ispwInCalTable][f]: |
538 - | matched, mymatch = sloppyMatch(t, timerangeListTimes, solutionTimeThresholdSeconds, |
539 - | myprint=debugSloppyMatch, whichone=True) |
540 - | if (matched): |
541 - | uTPFPS.append(t) |
563 + | found = 0 |
564 + | for timestamp in uniqueTimesPerFieldPerSpw[ispwInCalTable][f]: |
565 + | if timestamp in uniqueTimes: # 2024Aug28 |
566 + | myUniqueTimeIndex = uniqueTimes.index(timestamp) # 2024Aug28 |
567 + | matched, mymatch = sloppyMatch(timestamp, timerangeListTimes, solutionTimeThresholdSeconds, myUniqueTimeIndex, scansToPlot, # 2024Aug28 |
568 + | scansForUniqueTimes, myprint=debugSloppyMatch, whichone=True) |
569 + | # when there are multiple solutions per scan (as in SDSKY_PS tables), the value of mymatch will always be the first one, so cannot use it |
570 + | mymatch = myUniqueTimeIndex # 2024Aug28 |
571 + | else: |
572 + | # Not including the mytime, scansToPlot and scansForUniqueTimes in the following call will lead to the legend |
573 + | # showing all timestamps, not merely the subset that was selected by the (optional) scans parameter. That is |
574 + | # why we need to compare each myUniqueTimeIndex. You cannot simply use mytime as it will always be the final |
575 + | # time in the parent loop when you are in the overlay='time' scenario. But we fail over to this method. |
576 + | if debug: |
577 + | print("%.3f not in uniqueTimes, failing over to simpler method" % (timestamp)) |
578 + | matched, mymatch = sloppyMatch(timestamp, timerangeListTimes, solutionTimeThresholdSeconds, myprint=debugSloppyMatch, whichone=True) # 2024Aug27 regression |
579 + | if (matched and mymatch in timerangeList): |
580 + | uTPFPS.append(timestamp) |
542 581 | uTPFPStimerange.append(mymatch) |
543 - | |
582 + | expected = len(uniqueTimesPerFieldPerSpw[ispwInCalTable][f]) |
583 + | if found < expected: |
584 + | statement = "plotbandpass drawOverlayTimeLegends() found %d/%d time matches for spw%d field%d" % (found,expected,myspw,f) |
585 + | casalogPost(debug,statement) |
586 + | # sort the timeranges so the text labels will be increasing |
544 587 | idx = np.argsort(uTPFPS) |
545 588 | uTPFPStimerange = np.array(uTPFPStimerange)[idx] |
546 589 | uTPFPS = np.sort(uTPFPS) |
547 590 | timeFormat = 3 # HH:MM:SS |
548 591 | maxTimesAcross = maxTimesAcrossTheTop |
549 592 | if (firstFrame == 111): |
550 593 | maxTimesAcross -= 2 |
594 + | |
551 595 | for a in range(len(uTPFPS)): |
552 596 | legendString = utstring(uTPFPS[a],timeFormat) |
553 597 | if (debug): print("----> Defined legendString: %s" % (legendString)) |
554 598 | if (a==0): |
555 599 | pb.text(xstartTitle-0.03, ystartOverlayLegend, 'UT',color='k',fontsize=mysize, |
556 600 | transform=pb.gcf().transFigure) |
557 601 | if (a < maxTimesAcross): |
558 602 | x0 = xstartTitle + (a*timeHorizontalSpacing) |
559 603 | y0 = ystartOverlayLegend |
560 604 | else: |
561 605 | # start going down the righthand side |
562 606 | x0 = xstartTitle + (maxTimesAcross*timeHorizontalSpacing) |
563 607 | y0 = ystartOverlayLegend-(a-maxTimesAcross)*antennaVerticalSpacing |
564 - | # for tlt in timerangeListTimes: |
565 608 | if (True): |
566 609 | if (debug): |
567 610 | print("3)checking time %d, len(uTPFPS)=%d" % (a,len(uTPFPS))) |
568 611 | if (sloppyMatch(uTPFPS[a],timerangeListTimes, |
569 612 | solutionTimeThresholdSeconds, |
570 613 | mytime, scansToPlot, scansForUniqueTimes, |
571 614 | myprint=debugSloppyMatch)): |
572 615 | myUniqueTime = uTPFPS[a] |
573 616 | if (debug): |
574 617 | print("3)setting myUniqueTime to %d" % (myUniqueTime)) |
575 618 | |
576 619 | if (debug): print("----> Drawing legendString: %s" % (legendString)) |
577 620 | if ((len(fieldsToPlot) > 1 or len(timerangeList) > 1) and overlayAntennas==False): |
578 621 | # having overlayAntennas==False here will force all time labels to be black (as desired) |
579 622 | if (debug): |
623 + | print("len(uTPFPS)=%d, a=%d, len(myUniqueColor)=%d, overlayColors[%d]=%s" % (len(uTPFPS),a,len(myUniqueColor),timerangeList[allTimeranges.index(uTPFPStimerange[a])],str(overlayColors[timerangeList[allTimeranges.index(uTPFPStimerange[a])]]))) |
580 624 | print("len(uTPFPS)=%d, a=%d, len(myUniqueColor)=%d" % (len(uTPFPS),a,len(myUniqueColor))) |
581 - | # pb.text(x0, y0, legendString,color=overlayColors[timerangeList[a]],fontsize=mysize, |
582 - | # transform=pb.gcf().transFigure) |
583 - | if (debug): |
584 - | print("len(uTPFPStimerange)=%d, a=%d, len(myUniqueColor)=%d" % (len(uTPFPStimerange),a,len(myUniqueColor))) |
585 625 | |
586 - | uTPFPStimerangeValue = uTPFPStimerange[a] |
587 - | try: |
588 - | timerangeListIndex = allTimeranges.index(uTPFPStimerangeValue) |
589 - | except ValueError: |
590 - | casalogPost(debug, "uTPFPStimerangeValue = {} is not in the list allTimeranges".format(uTPFPStimerangeValue)) |
591 - | casalogPost(debug, "Setting timerangeListIndex = 0. This will change the overlay colors index to a default value. The time labels might not look as expected.") |
592 - | casalogPost(debug, "It is possible that certain antennas had a slightly different timestamp.") |
593 - | timerangeListIndex = 0 |
594 - | |
595 - | overlayColorsIndex = timerangeList[timerangeListIndex] |
596 - | pb.text(x0, y0, legendString, color=overlayColors[overlayColorsIndex], fontsize=mysize, transform=pb.gcf().transFigure) |
597 - | |
598 - | if (debug): |
599 - | print("done text") |
600 - | else: |
626 + | mycolor = overlayColors[uTPFPStimerange[a]] # color based on the subset of timeranges to be plotted |
627 + | pb.text(x0, y0, legendString, color=mycolor, fontsize=mysize, transform=pb.gcf().transFigure) |
628 + | else: # only 1 spectrum, or overlayAntennas==True means use all black labels |
601 629 | pb.text(x0, y0, legendString, fontsize=mysize, transform=pb.gcf().transFigure) |
602 630 | |
603 631 | def lineNumber(): |
604 632 | """Returns the current line number in our program.""" |
605 633 | return inspect.currentframe().f_back.f_lineno |
606 634 | |
607 635 | def drawAtmosphereAndFDM(showatm, showtsky, atmString, subplotRows, mysize, TebbSky, |
608 636 | TebbSkyImage,plotrange, xaxis, atmchan, atmfreq, transmission, |
609 637 | subplotCols, showatmPoints,xframe, channels,LO1,atmchanImage, |
610 638 | atmfreqImage,transmissionImage, firstFrame,showfdm,nChannels,tableFormat, |
705 733 | casalog.post(mystring) |
706 734 | if (debug): print(mystring) |
707 735 | |
708 736 | def computeHighestSpwIndexInSpwsToPlotThatHasCurrentScan(spwsToPlot, scansToPlotPerSpw, scan): |
709 737 | highestSpwIndex = -1 |
710 738 | for i,spw in enumerate(spwsToPlot): |
711 739 | if (scan in scansToPlotPerSpw[spw]): |
712 740 | highestSpwIndex = i |
713 741 | return(highestSpwIndex) |
714 742 | |
743 + | def madOfDiff(solution): |
744 + | """ |
745 + | This function is used to decide which of two curves has more scatter, and hence should |
746 + | be plotted first (i.e. shown in the background) when overlaying two solutions. |
747 + | Added as part of CAS-9474 to do a better job of the selection |
748 + | """ |
749 + | if (len(solution) < 4): |
750 + | return MAD(np.diff(solution)) |
751 + | else: |
752 + | start = int(len(solution)/4) |
753 + | stop = int(len(solution)*3/4) |
754 + | ydata = np.array(solution[start:stop+1]) |
755 + | return MAD(np.diff(ydata)) |
756 + | |
715 757 | def run_with_old_pyplot_style(func): |
716 758 | """ |
717 759 | CAS-12786: this decorator is introduced here to have a single entry point before/after |
718 760 | which the style sheet of matplotlib/pylab can be tuned and restored before returning. |
719 761 | |
720 762 | The motivation is that plotbandpass plots are heavily dependent on default image size, |
721 763 | grid style, ticks frequency, etc., as it used to be in matplotlib 1.1.0 (used in |
722 764 | CASA 5.x). In more recent matplotlib (3.1.1 is currently used in CASA 6.x) that style |
723 765 | can be reproduced using the style sheet 'classic'. This is the quickest and simplest way |
724 766 | to produce plotbandpass plots that look like those of CASA 5 |
1309 1351 | net_sideband = mytb.getcol('NET_SIDEBAND') |
1310 1352 | measFreqRef = mytb.getcol('MEAS_FREQ_REF') |
1311 1353 | originalSpw_casa33 = list(range(len(measFreqRef))) |
1312 1354 | chanFreqGHz_casa33 = [] # used by showFDM |
1313 1355 | for i in originalSpw_casa33: |
1314 1356 | # They array shapes can vary. |
1315 1357 | chanFreqGHz_casa33.append(1e-9 * mytb.getcell('CHAN_FREQ',i)) |
1316 1358 | mytb.close() |
1317 1359 | except: |
1318 1360 | print("2) Could not open the associated measurement set tables (%s). Will not translate antenna names." % (msName)) |
1361 + | cal_scans = None #### added 2024Aug |
1319 1362 | else: # 3.4 |
1320 1363 | tableFormat = 34 |
1321 1364 | cal_desc_id = mytb.getcol('SPECTRAL_WINDOW_ID') |
1322 1365 | cal_scans = mytb.getcol('SCAN_NUMBER') |
1323 1366 | unique_cal_scans = np.unique(cal_scans) |
1324 1367 | cal_scans_per_spw = {} |
1325 1368 | for myspw in np.unique(cal_desc_id): |
1326 1369 | cal_scans_per_spw[myspw] = np.unique(cal_scans[np.where(myspw == cal_desc_id)[0]]) |
1327 1370 | if (debug): |
1328 1371 | print("spw %d: scans %s" % (myspw,str(cal_scans_per_spw[myspw]))) |
1435 1478 | missingFrequencyWidth = originalChannelStart[ictr]*(mymsmd.chanfreqs(i)[-1]-mymsmd.chanfreqs(i)[0])/(nchan-1) |
1436 1479 | else: |
1437 1480 | missingFrequencyWidth = 0 |
1438 1481 | if (missingFrequencyWidth > 0): |
1439 1482 | if (DEBUG): |
1440 1483 | print("Correcting for channels flagged prior to running bandpass by %f GHz" % (missingFrequencyWidth*1e-9)) |
1441 1484 | newfreqs = 1e-9*(mymsmd.chanfreqs(i)) + missingFrequencyWidth*1e-9 |
1442 1485 | if debug: print("Appending onto chanFreqGHz: %s" % (str(newfreqs))) |
1443 1486 | chanFreqGHz.append(newfreqs) |
1444 1487 | |
1488 + | # the sort order of this variable is based on tb.getcol('SPECTRAL_WINDOW_ID') which is usually (always?) in increasing order |
1445 1489 | uniqueSpwsInCalTable = np.unique(cal_desc_id) |
1446 1490 | |
1447 1491 | # initial calculation for final message if not all spws appear with overlay='antenna' |
1448 1492 | uniqueTimes = sloppyUnique(np.unique(times), 1.0) |
1449 1493 | nUniqueTimes = len(uniqueTimes) |
1450 1494 | if (nUniqueTimes == 1): |
1451 1495 | solutionTimeSpread = 0 |
1452 1496 | else: |
1453 1497 | solutionTimeSpread = np.max(uniqueTimes)-np.min(uniqueTimes) |
1454 - | casalogPost(debug,"Found solutions with %d unique times (within a threshold of 1.0 second)." % (nUniqueTimes)) |
1455 - | |
1456 - | print("Median difference between solution times = %f sec" % (np.median(np.diff(times)))) |
1457 - | uniqueTimes = sloppyUnique(np.unique(times), solutionTimeThresholdSeconds) |
1498 + | casalogPost(debug,"Found solutions with %d unique times across all spws and fields (within a threshold of 1.0 second)." % (nUniqueTimes)) |
1499 + | |
1500 + | casalogPost(True,"Median difference between solution times = %f sec" % (np.median(np.diff(uniqueTimes)))) |
1501 + | if cal_scans is None or VisCal == 'SDSKY_PS': ### added 2024Aug |
1502 + | # this will show the spectrum for every solution / integration (appropriate for old Tsys tables and SDsky spectra) |
1503 + | uniqueTimes = sloppyUnique(np.unique(times), solutionTimeThresholdSeconds) ### indented 2024Aug |
1504 + | else: ### added 2024Aug |
1505 + | # this will show one spectrum per scan (appropriate for Tsys tables in CASA 3.4 onward) |
1506 + | uniqueTimes = sloppyUnique(times, solutionTimeThresholdSeconds, cal_scans) ### added 2024Aug |
1458 1507 | nUniqueTimes = len(uniqueTimes) |
1459 1508 | if (nUniqueTimes == 1): |
1460 1509 | casalogPost(debug,"Found solutions with %d unique time (within a threshold of %d seconds)." % (nUniqueTimes,solutionTimeThresholdSeconds)) |
1461 1510 | else: |
1462 1511 | casalogPost(debug,"Found solutions with %d unique times (within a threshold of %d seconds)." % (nUniqueTimes,solutionTimeThresholdSeconds)) |
1463 1512 | |
1464 1513 | scansForUniqueTimes = [] |
1465 1514 | if (tableFormat >= 34): |
1466 1515 | if (len(unique_cal_scans) == 1): |
1467 1516 | casalogPost(debug,"Found solutions with %d unique scan number %s" % (len(unique_cal_scans), str(unique_cal_scans))) |
1485 1534 | nFields = len(uniqueFields) |
1486 1535 | spwlist = [] |
1487 1536 | uniqueTimesPerFieldPerSpw = [] |
1488 1537 | for s in uniqueSpwsInCalTable: |
1489 1538 | uniqueTimesPerField = [] |
1490 1539 | for f in uniqueFields: |
1491 1540 | timelist = [] |
1492 1541 | for row in range(len(fields)): |
1493 1542 | if (fields[row] == f and cal_desc_id[row] == s): |
1494 1543 | if (sloppyMatch(times[row], timelist, solutionTimeThresholdSeconds) == False): |
1544 + | # if this time is not already in the list for this spw/field combination, then append it |
1495 1545 | timelist.append(times[row]) |
1496 1546 | spwlist.append(cal_desc_id) |
1497 1547 | uniqueTimesPerField.append(timelist) |
1498 1548 | uniqueTimesPerFieldPerSpw.append(uniqueTimesPerField) |
1499 1549 | |
1500 1550 | if (debug): print("about to call casalogPost") |
1501 1551 | casalogPost(debug,displayTimesArray([[uniqueTimes]])) |
1502 1552 | |
1503 1553 | # Parse the spws to plot from the command line |
1504 1554 | if (spw==''): |
1517 1567 | break |
1518 1568 | elif (token.find('~')>0): |
1519 1569 | (start,finish) = token.split('~') |
1520 1570 | spwsToPlot += list(range(int(start),int(finish)+1)) |
1521 1571 | else: |
1522 1572 | spwsToPlot.append(int(token)) |
1523 1573 | elif (type(spw) == list): |
1524 1574 | spwsToPlot = np.sort(spw) |
1525 1575 | else: |
1526 1576 | spwsToPlot = [spw] |
1577 + | # note that spwsToPlot will not necessarily be in increasing order, e.g. if the user specified them out-of-order |
1527 1578 | |
1528 1579 | casalogPost(debug,"%d spw%s in the solution = %s" % (len(uniqueSpwsInCalTable), plural(uniqueSpwsInCalTable), str(uniqueSpwsInCalTable))) |
1529 1580 | keepSpwsToPlot = spwsToPlot[:] |
1530 1581 | for myspw in spwsToPlot: |
1531 1582 | if (myspw not in uniqueSpwsInCalTable): |
1532 1583 | print("WARNING: spw %d is not in the solution. Removing it from the list to plot." % (myspw)) |
1533 1584 | print("Available spws = ", uniqueSpwsInCalTable) |
1534 1585 | keepSpwsToPlot.remove(myspw) |
1535 1586 | if (ctsys.compare_version('>=',[4,1,0]) and mymsmd != ''): |
1536 1587 | # # nonwvrspws = list(set(range(mymsmd.nspw())).difference(set(mymsmd.wvrspws()))) |
1886 1937 | if (myscan not in scansForUniqueTimes): |
1887 1938 | print(("No rows for scan %d, only " % (myscan), np.unique(scansForUniqueTimes))) |
1888 1939 | return |
1889 1940 | timerangeOffset = scansForUniqueTimes.index(myscan) |
1890 1941 | timerangeList = np.array(timerangeList) + timerangeOffset |
1891 1942 | if (debug): print(("Since both timeranges and scans was specified, generated new effective timerangeList: ", timerangeList)) |
1892 1943 | if (max(timerangeList) >= len(uniqueTimes)): |
1893 1944 | print("Invalid timerange. Solution has %d times (%d~%d)" % (len(uniqueTimes),0,len(uniqueTimes)-1)) |
1894 1945 | return |
1895 1946 | timerangeListTimes = np.array(uniqueTimes)[timerangeList] |
1947 + | if debug: |
1948 + | print("%d timerangeListTimes to be plotted" % (len(timerangeListTimes))) |
1896 1949 | timerangeListTimesString = mjdsecArrayToUTString(timerangeListTimes) |
1897 1950 | if (tableFormat == 33 or scansForUniqueTimes == []): |
1898 1951 | # SMA data with scan numbers of -1 has empty list for scansForUniqueTimes |
1899 1952 | scansToPlot = [] |
1900 1953 | if (scans != ''): |
1901 1954 | print("Selection by scan is not possible for this dataset.") |
1902 1955 | return |
1903 1956 | else: |
1904 - | if (debug): print("scansForUniqueTimes = %s" % (str(scansForUniqueTimes))) |
1957 + | if (debug): print("A)scansForUniqueTimes = %s" % (str(scansForUniqueTimes))) |
1905 1958 | scansToPlot = np.array(scansForUniqueTimes)[timerangeList] |
1906 1959 | if (np.unique(scansToPlot)[0] == -1): |
1907 1960 | # scan numbers are not correct in this new-style cal table |
1908 1961 | scansToPlot = [] |
1909 1962 | if (scans != ''): |
1910 1963 | print("Selection by scan number is not possible with this dataset.") |
1911 1964 | return |
1912 1965 | if (scans != '' and scans != []): |
1913 1966 | if (type(scans) == list): |
1914 1967 | scansToPlot = scans |
1915 1968 | elif (type(scans) == str): |
1916 1969 | scansToPlot = [int(a) for a in scans.split(',')] |
1917 1970 | else: |
1918 1971 | scansToPlot = [scans] |
1919 1972 | for scan in scansToPlot: |
1920 1973 | if (scan not in scansForUniqueTimes): |
1921 1974 | print("Scan %d is not in any solution" % (scan)) |
1922 1975 | return |
1976 + | if debug: print("scans to plot: %s" % (str(scansToPlot))) |
1923 1977 | scansToPlotPerSpw = {} |
1924 - | for myspw in np.unique(cal_desc_id): |
1978 + | casalogPost(debug,"originalSpwsToPlot: %s" % (originalSpwsToPlot)) |
1979 + | # for myspw in np.unique(cal_desc_id): ### removed 2024Aug |
1980 + | for myspw in originalSpwsToPlot: ### added 2024Aug |
1925 1981 | scansToPlotPerSpw[myspw] = [] |
1926 - | for scan in scansToPlot: |
1927 - | for myspw in np.unique(cal_desc_id): |
1928 - | if (scan in cal_scans_per_spw[myspw]): |
1929 - | scansToPlotPerSpw[myspw].append(scan) |
1982 + | if tableFormat == 34 and scansForUniqueTimes != []: ### added 2024Aug |
1983 + | scansToPlotRevised = [] ### added 2024Aug |
1984 + | timerangeListTimesRevised = [] ### added 2024Aug |
1985 + | for scan in scansToPlot: ### indented 2024Aug |
1986 + | # for myspw in np.unique(cal_desc_id): ### removed 2024Aug |
1987 + | for myspw in originalSpwsToPlot: ### added 2024Aug |
1988 + | if (scan in cal_scans_per_spw[myspw]): ### indented 2024Aug |
1989 + | scansToPlotPerSpw[myspw].append(scan) ### indented 2024Aug |
1990 + | scansToPlotRevised.append(scan) #### added 2024Aug |
1991 + | idx = np.where(cal_scans == scan) #### added 2024Aug |
1992 + | timerangeListTimesRevised += list(np.unique(times[idx])) #### added 2024Aug |
1993 + | scansToPlot = np.unique(scansToPlotRevised) #### added 2024Aug |
1994 + | # timerangeListTimes = np.unique(timerangeListTimesRevised) #### added 2024Aug |
1995 + | timerangeListTimes = np.array(uniqueTimes)[timerangeList] #### added 2024Aug after the Aug22 benchmark |
1996 + | casalogPost(debug, "%d timerangeListTimes (after filtering for spw) = %s" % (len(timerangeListTimes), timerangeListTimes)) #### add 2024Aug08 |
1997 + | casalogPost(debug, "scansToPlot (after filtering for spw) = %s" % (scansToPlot)) #### added 2024Aug |
1930 1998 | |
1931 1999 | # remove spws that do not have any scans to be plotted |
1932 2000 | # but only for tables that have a scan number column, and not filled with all -1 |
1933 2001 | if (tableFormat > 33 and scansForUniqueTimes != []): |
1934 2002 | for myspw in np.unique(cal_desc_id): |
1935 - | if (debug): |
1936 - | print("scans to plot for spw %d: %s" % (myspw, scansToPlotPerSpw[myspw])) |
1937 - | if (scansToPlotPerSpw[myspw] == []): |
1938 - | indexDelete = np.where(spwsToPlot==myspw)[0] |
1939 - | if (len(indexDelete) > 0): |
1940 - | spwsToPlot = np.delete(spwsToPlot, indexDelete[0]) |
1941 - | print("spws to plot = ", spwsToPlot) |
2003 + | if myspw in scansToPlotPerSpw: #### added 2024Aug |
2004 + | if (scansToPlotPerSpw[myspw] == []): #### indented 2024Aug |
2005 + | indexDelete = np.where(spwsToPlot==myspw)[0] #### indented 2024Aug |
2006 + | if (len(indexDelete) > 0): #### indented 2024Aug |
2007 + | spwsToPlot = np.delete(spwsToPlot, indexDelete[0]) #### indented 2024Aug |
2008 + | elif (debug): |
2009 + | print("scans to plot for spw %d: %s" % (myspw, scansToPlotPerSpw[myspw])) |
2010 + | print("spws to plot (sorted) = ", sorted(spwsToPlot)) |
1942 2011 | casalogPost(debug,"scans to plot: %s" % (str(scansToPlot))) |
1943 2012 | casalogPost(debug,"UT times to plot: %s" % (timerangeListTimesString)) |
1944 2013 | casalogPost(debug,"Corresponding time IDs (0-based): %s" % (str(timerangeList))) |
1945 2014 | if (len(timerangeListTimes) > len(np.unique(scansToPlot))): |
1946 2015 | # fix for CAS-9474 |
1947 2016 | uniqueScansToPlot, idx = np.unique(scansToPlot, return_index=True) |
1948 2017 | if (len(uniqueScansToPlot) < len(scansToPlot)): |
1949 2018 | # If the solution time for one spw differs by more than solutionTimeThresholdSeconds from |
1950 2019 | # another spw, then we will get 2 identical entries for the same scan, and thus duplicate |
1951 2020 | # plots. So, remove one. |
2193 2262 | if (showatmfield not in fieldsToPlot): |
2194 2263 | print("The showatmfield (%d=%s) is not in the list of fields to plot: %s" % (showatmfield, showatmfieldName, str(fieldsToPlot))) |
2195 2264 | return() |
2196 2265 | |
2197 2266 | for i in fieldsToPlot: |
2198 2267 | match = np.where(i==uniqueFields)[0] |
2199 2268 | if (len(match) < 1 and bOverlay): |
2200 2269 | match = np.where(i==uniqueFields2)[0] |
2201 2270 | fieldIndicesToPlot.append(match[0]) |
2202 2271 | |
2203 - | casalogPost(debug,"spws to plot = %s" % (str(spwsToPlot))) |
2272 + | casalogPost(debug,"spws to plot (sorted) = %s" % (str(sorted(spwsToPlot)))) |
2204 2273 | casalogPost(debug,"Field IDs to plot: %s" % (str(fieldsToPlot))) |
2205 2274 | |
2206 2275 | redisplay = False |
2207 2276 | myap = 0 # this variable is necessary to make the 'b' option work for |
2208 2277 | # subplot=11, yaxis=both. It keeps track of whether 'amp' or |
2209 2278 | # 'phase' was the first plot on the page. |
2210 2279 | |
2211 2280 | # I added pb.ion() because Remy suggested it. |
2212 2281 | if (interactive): |
2213 2282 | pb.ion() # This will open a new window if not present. |
2237 2306 | if (corr_type_string == []): |
2238 2307 | return() |
2239 2308 | polsToPlot = checkPolsToPlot(polsToPlot, corr_type_string, debug) |
2240 2309 | if (polsToPlot == []): |
2241 2310 | return() |
2242 2311 | # Here we are only plotting one BPOLY solution, no overlays implemented. |
2243 2312 | overlayAntennas = False |
2244 2313 | # rows in the table are: antennas 0..nAnt for first spw, antennas 0..nAnt |
2245 2314 | # for 2nd spw... |
2246 2315 | pagectr = 0 |
2316 + | if debug: |
2317 + | print("Setting pages to blank list") |
2247 2318 | pages = [] |
2248 2319 | xctr = 0 |
2249 2320 | newpage = 1 |
2250 2321 | while (xctr < len(antennasToPlot)): |
2251 2322 | xant = antennasToPlot[xctr] |
2252 2323 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
2253 2324 | spwctr = 0 |
2254 2325 | spwctrFirstToPlot = spwctr |
2255 2326 | while (spwctr < len(spwsToPlot)): |
2256 2327 | ispw = spwsToPlot[spwctr] |
2257 2328 | mytime = 0 |
2258 2329 | while (mytime < nUniqueTimes): |
2259 2330 | if (len(uniqueTimes) > 0 and (mytime not in timerangeList)): |
2260 2331 | if (debug): |
2261 2332 | print("@@@@@@@@@@@@@@@ Skipping mytime=%d" % (mytime)) |
2262 2333 | mytime += 1 |
2263 2334 | continue |
2264 2335 | if (newpage == 1): |
2265 2336 | pages.append([xctr,spwctr,mytime,0]) |
2266 - | # print("appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,0)) |
2337 + | if debug: |
2338 + | print("top: appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,0)) |
2267 2339 | newpage = 0 |
2268 2340 | antennaString = 'Ant%2d: %s, ' % (xant,antstring) |
2269 2341 | for index in range(nRows): |
2270 2342 | # Find this antenna, spw, and timerange combination in the table |
2343 + | if tableFormat >= 34: ### added 2024Aug |
2344 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
2345 + | else: ### added 2024Aug |
2346 + | scansToPlotHere = scansToPlot ### added 2024Aug |
2271 2347 | if (xant==ant[index] and sloppyMatch(uniqueTimes[mytime],times[index],solutionTimeThresholdSeconds, |
2272 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, |
2348 + | mytime, scansToPlotHere, scansForUniqueTimes, ### modified 2024Aug |
2273 2349 | myprint=debugSloppyMatch) and |
2274 2350 | (ispw == cal_desc_id[index]) and (fields[index] in fieldsToPlot)): |
2275 2351 | fieldIndex = np.where(fields[index] == uniqueFields)[0] |
2276 2352 | if (type(fieldIndex) == list or type(fieldIndex) == np.ndarray): |
2277 2353 | fieldIndex = fieldIndex[0] |
2278 2354 | validDomain = [frequencyLimits[0,index], frequencyLimits[1,index]] |
2279 2355 | if (msFound): |
2280 2356 | fieldString = msFields[uniqueFields[fieldIndex]] |
2281 2357 | else: |
2282 2358 | fieldString = str(field) |
2283 2359 | timeString = ', t%d/%d %s' % (mytime,nUniqueTimes-1,utstring(uniqueTimes[mytime],3)) |
2284 2360 | if (scansForUniqueTimes != []): |
2285 2361 | if (scansForUniqueTimes[mytime]>=0): |
2286 2362 | timeString = ', scan%d %s' % (scansForUniqueTimes[mytime],utstring(uniqueTimes[mytime],3)) |
2287 2363 | if ((yaxis.find('amp')>=0 or amplitudeWithPhase) and myap==0): |
2288 2364 | xframe += 1 |
2289 2365 | myUniqueColor = [] |
2290 2366 | if (debug): |
2291 2367 | print("v) incrementing xframe to %d" % xframe) |
2292 - | adesc = pb.safe_pb_subplot(xframe) |
2368 + | adesc = safe_pb_subplot(xframe) |
2293 2369 | previousSubplot = xframe |
2294 2370 | if (ispw==originalSpw[ispw]): |
2295 2371 | # all this was added mistakenly here. If it causes a bug, remove it. |
2296 2372 | if (overlayTimes and len(fieldsToPlot) > 1): |
2297 2373 | indices = fstring = '' |
2298 2374 | for f in fieldIndicesToPlot: |
2299 2375 | if (f != fieldIndicesToPlot[0]): |
2300 2376 | indices += ',' |
2301 2377 | fstring += ',' |
2302 2378 | indices += str(uniqueFields[f]) |
2341 2417 | SetNewXLimits([plotrange[0],plotrange[1]]) |
2342 2418 | if (plotrange[2] != 0 or plotrange[3] != 0): |
2343 2419 | SetNewYLimits([plotrange[2],plotrange[3]]) |
2344 2420 | xlim=pb.xlim() |
2345 2421 | ylim=pb.ylim() |
2346 2422 | ResizeFontsSetGrid(adesc,mysize) |
2347 2423 | if (yaxis.lower().find('db')>=0): |
2348 2424 | pb.ylabel('Amplitude (dB)', size=mysize) |
2349 2425 | else: |
2350 2426 | pb.ylabel('Amplitude', size=mysize) |
2351 - | pb.xlabel('Frequency (GHz)', size=mysize) |
2427 + | # pb.xlabel('Frequency (GHz)', size=mysize) |
2428 + | pb.xlabel('Frequency (GHz) (%d channels)'%(len(frequenciesGHz[index])), size=mysize) # 2024Aug |
2352 2429 | if (xframe == firstFrame): |
2353 2430 | DrawBottomLegendPageCoords(msName, uniqueTimes[mytime], mysize, figfile) |
2354 2431 | pb.text(xstartTitle, ystartTitle, |
2355 2432 | '%s (degamp=%d, degphase=%d)'%(caltableTitle,nPolyAmp[index]-1, |
2356 2433 | nPolyPhase[index]-1),size=mysize, |
2357 2434 | transform=pb.gcf().transFigure) |
2358 2435 | # draw polarization labels |
2359 2436 | x0 = xstartPolLabel |
2360 2437 | y0 = ystartPolLabel |
2361 2438 | for p in range(nPolarizations): |
2362 2439 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
2363 2440 | pb.text(x0, y0-0.03*subplotRows*p, corrTypeToString(corr_type[p])+'', |
2364 2441 | color=pcolor[p],size=mysize, transform=pb.gca().transAxes) |
2365 2442 | if (xframe == 111 and amplitudeWithPhase): |
2366 2443 | if (len(figfile) > 0): |
2367 2444 | # We need to make a new figure page |
2368 2445 | plotfiles.append(makeplot(figfile,msFound,msAnt, |
2369 - | overlayAntennas,pages,pagectr, |
2370 - | density,interactive,antennasToPlot, |
2371 - | spwsToPlot,overlayTimes,overlayBasebands, |
2446 + | overlayAntennas,pages,pagectr, |
2447 + | density,interactive,antennasToPlot, |
2448 + | spwsToPlot,overlayTimes,overlayBasebands, |
2372 2449 | 0,xant,ispw,subplot,resample, |
2373 - | debug,figfileSequential,figfileNumber)) |
2450 + | debug,figfileSequential,figfileNumber)) |
2374 2451 | figfileNumber += 1 |
2375 2452 | donetime = time.time() |
2376 2453 | if (interactive): |
2377 2454 | pb.draw() |
2378 2455 | # # myinput = raw_input("(%.1f sec) Press return for next page (b for backwards, q to quit): "%(donetime-mytimestamp)) |
2379 2456 | myinput = input("Press return for next page (b for backwards, q to quit): ") |
2380 2457 | else: |
2381 2458 | myinput = '' |
2382 2459 | skippingSpwMessageSent = 0 |
2383 2460 | mytimestamp = time.time() |
2393 2470 | mytime = pages[pagectr][PAGE_TIME] |
2394 2471 | myap = pages[pagectr][PAGE_AP] |
2395 2472 | xant = antennasToPlot[xctr] |
2396 2473 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
2397 2474 | ispw = spwsToPlot[spwctr] |
2398 2475 | redisplay = True |
2399 2476 | else: |
2400 2477 | pagectr += 1 |
2401 2478 | if (pagectr >= len(pages)): |
2402 2479 | pages.append([xctr,spwctr,mytime,1]) |
2403 - | # print("appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,1)) |
2480 + | print("appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,1)) |
2404 2481 | newpage = 0 |
2405 2482 | safe_pb_clf() |
2406 2483 | |
2407 2484 | if (yaxis.find('phase')>=0 or amplitudeWithPhase): |
2408 2485 | xframe += 1 |
2409 2486 | myUniqueColor = [] |
2410 2487 | # # print("w) incrementing xframe to %d" % xframe) |
2411 - | adesc = pb.safe_pb_subplot(xframe) |
2488 + | adesc = safe_pb_subplot(xframe) |
2412 2489 | previousSubplot = xframe |
2413 2490 | if (ispw==originalSpw[ispw]): |
2414 2491 | pb.title("%sspw%2d, field %d: %s%s" % (antennaString,ispw, |
2415 2492 | uniqueFields[fieldIndex],fieldString,timeString), size=titlesize) |
2416 2493 | else: |
2417 2494 | pb.title("%sspw%2d (%d), field %d: %s%s" % (antennaString,ispw,originalSpw[ispw], |
2418 2495 | uniqueFields[fieldIndex],fieldString,timeString), size=titlesize) |
2419 2496 | phaseSolutionX = calcChebyshev(polynomialPhase[index][0:nPolyPhase[index]], validDomain, frequenciesGHz[index]*1e+9) * 180/math.pi |
2420 2497 | phaseSolutionY = calcChebyshev(polynomialPhase[index][nPolyPhase[index]:2*nPolyPhase[index]], validDomain, frequenciesGHz[index]*1e+9) * 180/math.pi |
2421 2498 | if (nPolarizations == 1): |
2422 2499 | pb.plot(frequenciesGHz[index], phaseSolutionX, '%s%s'%(xcolor,bpolymarkstyle),markeredgewidth=markeredgewidth) |
2423 2500 | else: |
2424 2501 | pb.plot(frequenciesGHz[index], phaseSolutionX, '%s%s'%(xcolor,bpolymarkstyle), frequenciesGHz[index], phaseSolutionY, '%s%s'%(ycolor,bpolymarkstyle),markeredgewidth=markeredgewidth) |
2425 2502 | ResizeFontsSetGrid(adesc,mysize) |
2426 2503 | pb.ylabel('Phase (deg)', size=mysize) |
2427 - | pb.xlabel('Frequency (GHz)', size=mysize) |
2504 + | # pb.xlabel('Frequency (GHz)', size=mysize) |
2505 + | pb.xlabel('Frequency (GHz) (%d channels)'%(len(frequenciesGHz[index])), size=mysize) ### 2024Aug |
2428 2506 | if (plotrange[0] != 0 or plotrange[1] != 0): |
2429 2507 | SetNewXLimits([plotrange[0],plotrange[1]]) |
2430 2508 | if (plotrange[2] != 0 or plotrange[3] != 0): |
2431 2509 | SetNewYLimits([plotrange[2],plotrange[3]]) |
2432 2510 | if (amplitudeWithPhase and phase != ''): |
2433 2511 | if (phase[0] != 0 or phase[1] != 0): |
2434 2512 | SetNewYLimits(phase) |
2435 2513 | if (xframe == firstFrame): |
2436 2514 | pb.text(xstartTitle, ystartTitle, |
2437 2515 | '%s (degamp=%d, degphase=%d)'%(caltable, |
2444 2522 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
2445 2523 | pb.text(x0, y0-0.03*p*subplotRows, corrTypeToString(corr_type[p])+'', |
2446 2524 | color=pcolor[p],size=mysize, transform=pb.gca().transAxes) |
2447 2525 | |
2448 2526 | # end of 'for' loop over rows |
2449 2527 | redisplay = False |
2450 2528 | pb.subplots_adjust(hspace=myhspace, wspace=mywspace) |
2451 2529 | if (xframe == lastFrame): |
2452 2530 | if (len(figfile) > 0): |
2453 2531 | plotfiles.append(makeplot(figfile,msFound,msAnt, |
2454 - | overlayAntennas,pages,pagectr, |
2455 - | density,interactive,antennasToPlot, |
2456 - | spwsToPlot,overlayTimes,overlayBasebands, |
2532 + | overlayAntennas,pages,pagectr, |
2533 + | density,interactive,antennasToPlot, |
2534 + | spwsToPlot,overlayTimes,overlayBasebands, |
2457 2535 | 1,xant,ispw, |
2458 2536 | subplot,resample,debug, |
2459 - | figfileSequential,figfileNumber)) |
2537 + | figfileSequential,figfileNumber)) |
2460 2538 | figfileNumber += 1 |
2461 2539 | donetime = time.time() |
2462 2540 | if (interactive): |
2463 2541 | pb.draw() |
2464 2542 | # # myinput = raw_input("(%.1f sec) Press return for next page (b for backwards, q to quit): "%(donetime-mytimestamp)) |
2465 2543 | myinput = input("Press return for next page (b for backwards, q to quit): ") |
2466 2544 | else: |
2467 2545 | myinput = '' |
2468 2546 | skippingSpwMessageSent = 0 |
2469 2547 | mytimestamp = time.time() |
2718 2796 | |
2719 2797 | madsigma = channeldiff # for option channeldiff>0, sets threshold for finding outliers |
2720 2798 | ampMin = LARGE_POSITIVE |
2721 2799 | ampMax = LARGE_NEGATIVE |
2722 2800 | PHASE_ABS_SUM_THRESHOLD = 2e-3 # in degrees, used to avoid printing MAD statistics for refant |
2723 2801 | |
2724 2802 | TDMisSecond = False |
2725 2803 | pagectr = 0 |
2726 2804 | drewAtmosphere = False |
2727 2805 | newpage = 1 |
2806 + | if debug: |
2807 + | print("Setting pages to blank list") |
2728 2808 | pages = [] |
2729 2809 | xctr = 0 |
2730 2810 | myap = 0 # determines whether an amp or phase plot starts the page (in the case of 'both') |
2731 2811 | # zero means amplitude, 1 means phase |
2732 2812 | redisplay = False |
2733 2813 | matchctr = 0 |
2734 2814 | myUniqueColor = [] |
2735 2815 | # for the overlay=antenna case, start by assuming the first antenna is not flagged |
2736 2816 | firstUnflaggedAntennaToPlot = 0 |
2737 2817 | lastUnflaggedAntennaToPlot = len(antennasToPlot) |
2827 2907 | if (overlaySpws or overlayBasebands): |
2828 2908 | groupByBaseband = True |
2829 2909 | if (groupByBaseband and overlaySpws==False and overlayBasebands==False): |
2830 2910 | showBasebandNumber = True |
2831 2911 | # Basic nested 'while' loop structure is: |
2832 2912 | # - antennas |
2833 2913 | # - baseband (if necessary) |
2834 2914 | # - spw |
2835 2915 | # - time |
2836 2916 | # - for i in rows |
2917 + | maxChannels = {}; maxChannels2 = {} |
2837 2918 | while (xctr < len(antennasToPlot)): |
2838 2919 | if (debug): print("at top of xctr loop: %d" % (xctr)) |
2839 2920 | xant = antennasToPlot[xctr] |
2840 2921 | bbctr = 0 |
2841 2922 | spwctr = 0 |
2842 2923 | spwctrFirstToPlot = 0 |
2843 2924 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
2844 2925 | alreadyPlottedAmp = False # needed for (overlay='baseband', yaxis='both') CAS-6477 |
2845 2926 | finalSpwWasFlagged = False # inserted on 22-Apr-2014 for g25.27 |
2846 2927 | if debug: |
2917 2998 | nUniqueTimes = len(uniqueTimes) |
2918 2999 | # currentSpwctr = spwctr # commented out on 2014-04-04 to match task for task01 regression |
2919 3000 | if (overlaySpws or overlayBasebands): |
2920 3001 | if (xctr >= firstUnflaggedAntennaToPlot): |
2921 3002 | if (debug): |
2922 3003 | print("xctr=%d >= firstUnflaggedAntennaToPlot=%d, decrementing spwctr to %d" % (xctr, firstUnflaggedAntennaToPlot,spwctr-1)) |
2923 3004 | spwctr -= 1 |
2924 3005 | |
2925 3006 | firstTimeMatch = -1 |
2926 3007 | finalTimeMatch = -1 # for CAS-7820 |
2927 - | while (mytime < nUniqueTimes): # start of enrmous 2470-line while loop |
3008 + | while (mytime < nUniqueTimes): # start of enormous 2470-line while loop |
2928 3009 | finalTimerangeFlagged = False # 04-Aug-2014 |
2929 3010 | if (debug): |
2930 - | print("mytime = %d < %d, uniqueTimes[mytime] = %s" % (mytime,nUniqueTimes,str(uniqueTimes[mytime]))) |
3011 + | print("at top of mytime loop: mytime = %d < %d" % (mytime,nUniqueTimes)) |
3012 + | print("%d timerangeListTimes" % (len(timerangeListTimes))) |
2931 3013 | print("timerangeList = %s" % (str(timerangeList))) |
2932 3014 | print("timerangeListTimes = %s" % (str(timerangeListTimes))) |
2933 3015 | print("debugSloppyMatch = %s" % (str(debugSloppyMatch))) |
2934 3016 | print(("solutionTimeThresholdSeconds = %s" % (str(solutionTimeThresholdSeconds)))) |
2935 3017 | # if ((scansToPlot == scansToPlotPerSpw[ispw]).all() == False and False): |
2936 3018 | # print " scansToPlot = ", scansToPlot |
2937 3019 | # print "scansToPlotPerSpw[%2d] = " % (ispw), scansToPlotPerSpw[ispw] |
2938 3020 | if (len(timerangeList) > 0 and |
2939 3021 | (sloppyMatch(uniqueTimes[mytime],timerangeListTimes,solutionTimeThresholdSeconds, |
2940 3022 | mytime, scansToPlot, scansForUniqueTimes, myprint=debugSloppyMatch)==False)): # task version |
2982 3064 | antennaString = 'Ant%2d: %s, ' % (xant,antstring) |
2983 3065 | if (overlayBasebands): |
2984 3066 | # Added on 7/29/2014 to fix infinite loop in uid___A002_X652932_X20fb bandpass |
2985 3067 | if (mytime == nUniqueTimes): |
2986 3068 | spwctr = len(spwsToPlot) |
2987 3069 | break |
2988 3070 | ispw = spwsToPlot[spwctr] |
2989 3071 | ispwInCalTable = list(uniqueSpwsInCalTable).index(ispw) |
2990 3072 | if (debug): |
2991 3073 | print("----------------------------- spwctr=%d, ispw set to %d, xctr=%d" % (spwctr,ispw,xctr)) |
2992 - | |
2993 - | # This used to be above the previous if/else block |
3074 + | # endif overlaySpws or overlayBasebands |
2994 3075 | if (newpage==1): |
2995 3076 | # add the current page (being created here) to the list |
3077 + | if (debug): |
3078 + | print("pages = ", pages) |
2996 3079 | pages.append([xctr,spwctr,mytime,0]) |
2997 3080 | if (debug): |
2998 - | print("top: appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,0)) |
3081 + | print("next: appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,0)) |
2999 3082 | newpage = 0 |
3083 + | if (ispw not in uniqueSpwsInCalTable): ##### added 2024Aug |
3084 + | print("spw %d is not in caltable=%s" % (ispw,uniqueSpwsInCalTable)) ##### added 2024Aug |
3085 + | return ##### added 2024Aug |
3086 + | if tableFormat > 33: ##### added 2024Aug |
3087 + | if ispw not in list(scansToPlotPerSpw.keys()): ##### added 2024Aug |
3088 + | print("ispw=%d not in %s" % (ispw,scansToPlotPerSpw)) ##### added 2024Aug |
3089 + | break ##### added 2024Aug |
3000 3090 | gplotx = [] |
3001 3091 | gploty = [] |
3002 3092 | channels = [] |
3003 3093 | xchannels = [] |
3004 3094 | ychannels = [] |
3005 3095 | frequencies = [] |
3006 3096 | xfrequencies = [] |
3007 3097 | yfrequencies = [] |
3008 3098 | channels2 = [] |
3009 3099 | xchannels2 = [] |
3026 3116 | if (overlayTimes or overlayAntennas or len(fieldsToPlot)>1 or |
3027 3117 | (nFields>1 and len(fieldlist)<nFields)): |
3028 3118 | # When there are multiple fields, then matching by scan causes the first |
3029 3119 | # matching solution to be displayed every time. So use the original method |
3030 3120 | # of matching by time until I think of something better. |
3031 3121 | sm = sloppyMatch(uniqueTimes[mytime],times[i],solutionTimeThresholdSeconds,myprint=False) |
3032 3122 | else: |
3033 3123 | if (overlayBasebands): |
3034 3124 | sTP = scansToPlot |
3035 3125 | else: |
3036 - | sTP = scansToPlotPerSpw[ispw] |
3126 + | if tableFormat > 33: ### added 2024Aug |
3127 + | sTP = scansToPlotPerSpw[ispw] ### indented 2024Aug |
3128 + | else: ### added 2024Aug |
3129 + | sTP = [] ### added 2024Aug |
3037 3130 | sm = sloppyMatch(uniqueTimes[mytime],times[i],solutionTimeThresholdSeconds, |
3038 3131 | mytime, sTP, scansForUniqueTimes, myprint=False) # au version |
3039 3132 | if ((ant[i]==xant) and (cal_desc_id[i]==ispw) and sm |
3040 3133 | and (mytime in timerangeList) # this test was added to support multiFieldInTimeOverlay |
3041 3134 | ): |
3042 3135 | if (debug): print("len(chanFreqGHz)=%d, ispw=%d" % (len(chanFreqGHz),ispw)) |
3043 3136 | if (msFound or tableFormat==34): |
3044 3137 | if (len(chanFreqGHz[ispw]) == 1): |
3045 3138 | if ((skippingSpwMessageSent & (1<<ispw)) == 0): |
3046 3139 | casalogPost(debug,"Skipping spw=%d because it has only 1 channel." % (ispw)) |
3079 3172 | gplotx.append(ggx[i][j]) |
3080 3173 | xchannels.append(j) |
3081 3174 | if (msFound or tableFormat==34): |
3082 3175 | xfrequencies.append(chanFreqGHz[ispw][j]) |
3083 3176 | if (nPolarizations == 2): |
3084 3177 | if (showflagged or (showflagged == False and flags[i][1][j]==0)): |
3085 3178 | gploty.append(ggy[i][j]) |
3086 3179 | ychannels.append(j) |
3087 3180 | if (msFound or tableFormat==34): |
3088 3181 | yfrequencies.append(chanFreqGHz[ispw][j]) |
3089 - | # end 'for i' |
3182 + | # end 'for i' over rows |
3183 + | |
3090 3184 | # if (not matchFound and newpage==0 and firstTimeMatch==-1): |
3091 3185 | if (not matchFound and newpage==0): |
3092 - | if (subplot==11 or (subplot!=11 and firstSpwMatch==-1 and firstTimeMatch==-1)): |
3186 + | # the first test below is the first of two fixes for CAS-13568 (prevent crash) |
3187 + | if (len(pages) > 1) and (subplot==11 or (subplot!=11 and firstSpwMatch==-1 and firstTimeMatch==-1)): |
3093 3188 | # Fix for CAS-7753 |
3094 3189 | # the firstTimeMatch part was needed for regression 65: different antennas having different solution times |
3095 3190 | newpage = 1 |
3191 + | if debug: print("setting pages to length=%d" % (len(pages)-1)) |
3096 3192 | pages = pages[:len(pages)-1] |
3097 3193 | myspw = originalSpw[ispw] |
3194 + | if myspw not in list(maxChannels.keys()): |
3195 + | maxChannels[myspw] = 0 # keep track to set x-axis label correctly |
3196 + | maxChannels2[myspw] = 0 # keep track to set x-axis label correctly |
3197 + | |
3198 + | if len(xchannels) > maxChannels[myspw]: |
3199 + | maxChannels[myspw] = len(xchannels) # keep track to set x-axis label correctly |
3200 + | if len(xchannels2) > maxChannels[myspw]: |
3201 + | maxChannels2[myspw] = len(xchannels2) # keep track to set x-axis label correctly |
3098 3202 | if (msFound): |
3099 3203 | if debug: |
3204 + | print("A) xchannels = ", xchannels) |
3100 3205 | print("myspw=%s" % (str(myspw))) |
3101 3206 | print("len(refFreq)=%d" % (len(refFreq))) |
3102 3207 | if (myspw >= len(refFreq)): |
3103 3208 | myspw = ispw |
3104 3209 | if (msFound and refFreq[myspw]*1e-9 > 60): |
3105 - | # Then this cannot be EVLA data. But I should really check the telescope name! |
3210 + | # Then this cannot be Band 1 or EVLA data. TODO: But I should really check the telescope name! |
3106 3211 | # if (refFreq[myspw]*1e-9 > np.mean(frequencies)): |
3107 - | if (refFreq[myspw]*1e-9 > np.mean(chanFreqGHz[ispw])): # this is safer (since frequencies might be []) |
3212 + | if (refFreq[myspw]*1e-9 > np.mean(chanFreqGHz[ispw])): # this is safer (since frequencies might be an empty list) |
3108 3213 | sideband = -1 |
3109 - | xlabelString = "%s LSB Frequency (GHz)" % refTypeToString(measFreqRef[myspw]) |
3214 + | # xlabelString = "%s LSB Frequency (GHz)" % refTypeToString(measFreqRef[myspw]) |
3215 + | xlabelString = "%s LSB Frequency (GHz) (%d channels)" % (refTypeToString(measFreqRef[myspw]),maxChannels[myspw]) ### 2024Aug |
3110 3216 | else: |
3111 3217 | sideband = +1 |
3112 - | xlabelString = "%s USB Frequency (GHz)" % refTypeToString(measFreqRef[myspw]) |
3218 + | # xlabelString = "%s USB Frequency (GHz)" % refTypeToString(measFreqRef[myspw]) |
3219 + | xlabelString = "%s USB Frequency (GHz) (%d channels)" % (refTypeToString(measFreqRef[myspw]),maxChannels[myspw]) ### 2024Aug |
3113 3220 | else: |
3114 3221 | sideband = -1 |
3115 - | xlabelString = "Frequency (GHz)" |
3222 + | # xlabelString = "Frequency (GHz)" |
3223 + | xlabelString = "Frequency (GHz) (%d channels)" % (maxChannels[myspw]) ### 2024Aug |
3116 3224 | if ((len(frequencies)>0) and (chanrange[1] > len(frequencies))): |
3117 3225 | print("Invalid chanrange (%d-%d) for spw%d in caltable1. Valid range = 0-%d" % (chanrange[0],chanrange[1],ispw,len(frequencies)-1)) |
3118 3226 | return() |
3119 3227 | pchannels = [xchannels,ychannels] |
3120 3228 | pfrequencies = [xfrequencies,yfrequencies] |
3121 3229 | gplot = [gplotx,gploty] |
3122 3230 | # We only need to compute the atmospheric transmission if: |
3123 3231 | # * we have been asked to show it, |
3124 3232 | # * there is a non-trivial number of channels, |
3125 3233 | # * the current field is the one for which we should calculate it (if times are being overlaied) |
3212 3320 | # difference between the requested Atm freq and actual Atm calcuation CAS-7715. No longer necessary after CAS-10228. |
3213 3321 | #chanDifference = atmfreq[0]-frequencies[0] |
3214 3322 | #chanImageDifference = atmfreqImage[0]-frequenciesImage[-1] |
3215 3323 | #print("signal SB difference = %f, image SB difference = %f" % (chanDifference, chanImageDifference)) |
3216 3324 | if showtsys: |
3217 3325 | TebbSkyImage = TsysImage |
3218 3326 | atmfreqImage = list(2*LO1 - np.array(atmfreqImage)) # + 2*chanDifference + chanImageDifference) # CAS-7715 adds final 2 terms |
3219 3327 | atmchanImage.reverse() |
3220 3328 | |
3221 3329 | if (overlayTimes): |
3222 - | atmString = 'PWV %.2fmm, airmass %.2f (field %d)' % (pwvmean,atmairmass,showatmfield) |
3330 + | atmString = 'PWV %.2fmm, airmass %.2f, maxAlt %.0fkm (field %d)' % (pwvmean,atmairmass,maxAltitude,showatmfield) |
3223 3331 | else: |
3224 - | atmString = 'PWV %.2fmm, airmass %.3f' % (pwvmean,atmairmass) |
3332 + | atmString = 'PWV %.2fmm, airmass %.3f, maxAlt %.0fkm' % (pwvmean,atmairmass,maxAltitude) |
3225 3333 | if (bOverlay): |
3226 3334 | for i in range(nRows2): |
3227 3335 | if (overlayTimes or overlayAntennas or len(fieldsToPlot)>1 or |
3228 3336 | (nFields>1 and len(fieldlist)<nFields)): |
3229 3337 | # Not having this path causes Tsys table overlays to behave like overlay='antenna,time' |
3230 3338 | # for caltable2. |
3231 3339 | sm = sloppyMatch(uniqueTimes2[mytime],times2[i],solutionTimeThresholdSeconds,myprint=False) |
3232 3340 | else: |
3233 3341 | if (mytime >= len(uniqueTimes2)): |
3234 3342 | # Fix for CAS-9474: avoid calling sloppyMatch because it will crash. |
3235 3343 | # Setting sm=False will result in an abort: "no amp data found in second solution." |
3236 3344 | sm = False |
3237 3345 | else: |
3346 + | if tableFormat >= 34: ### added 2024Aug |
3347 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
3348 + | else: ### added 2024Aug |
3349 + | scansToPlotHere = scansToPlot ### added 2024Aug |
3238 3350 | sm = sloppyMatch(uniqueTimes2[mytime],times2[i],solutionTimeThresholdSeconds, |
3239 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
3351 + | mytime, scansToPlotHere, scansForUniqueTimes, # au version ### modified 2024Aug |
3240 3352 | myprint=debugSloppyMatch) |
3241 3353 | if ((ant2[i]==xant) and (cal_desc_id2[i]==ispw) and sm |
3242 3354 | and (mytime in timerangeList) # added to match first caltable logic on 2014-04-09 |
3243 3355 | ): |
3244 3356 | if (fields2[i] in fieldsToPlot): |
3245 3357 | xflag2.append(flags2[i][0][:]) |
3246 3358 | yflag2.append(flags2[i][1][:]) |
3247 3359 | # With solint='2ch' or more, the following loop should not be over |
3248 3360 | # chanFreqGHz2 but over the channels in the solution. |
3249 3361 | for j in range(len(chanFreqGHz2[ispw])): |
3256 3368 | xfrequencies2.append(chanFreqGHz2[ispw][j]) |
3257 3369 | if (nPolarizations2 == 2): |
3258 3370 | if (showflagged or (showflagged == False and flags2[i][1][j]==0)): |
3259 3371 | gploty2.append(ggy2[i][j]) |
3260 3372 | ychannels2.append(j) |
3261 3373 | yfrequencies2.append(chanFreqGHz2[ispw][j]) |
3262 3374 | # end 'for i' |
3263 3375 | pchannels2 = [xchannels2,ychannels2] |
3264 3376 | pfrequencies2 = [xfrequencies2,yfrequencies2] |
3265 3377 | gplot2 = [gplotx2,gploty2] |
3378 + | # Need to rewrite the xlabel to show the total channel numbers from both caltables. Note that xaxis must be 'freq' for bOverlay |
3379 + | if (msFound and refFreq[myspw]*1e-9 > 60): |
3380 + | # Then this cannot be Band 1 or EVLA data. But I should really check the telescope name! |
3381 + | if (refFreq[myspw]*1e-9 > np.mean(chanFreqGHz2[ispw])): # this is safer (since frequencies might be []) |
3382 + | sideband = -1 |
3383 + | xlabelString = "%s LSB Frequency (GHz) (%d, %d channels)" % (refTypeToString(measFreqRef[myspw]),maxChannels[myspw],maxChannels2[myspw]) |
3384 + | else: |
3385 + | sideband = +1 |
3386 + | xlabelString = "%s USB Frequency (GHz) (%d, %d channels)" % (refTypeToString(measFreqRef[myspw]),maxChannels[myspw],maxChannels2[myspw]) |
3387 + | else: |
3388 + | sideband = -1 |
3389 + | xlabelString = "Frequency (GHz) (%d, %d channels)" % (maxChannels[myspw],maxChannels2[myspw]) |
3390 + | # endif bOverlay |
3266 3391 | |
3267 3392 | if (matchFound==False): |
3268 3393 | if ((overlayAntennas==False and overlaySpws==False and overlayBasebands==False) or |
3269 3394 | (overlayAntennas and xctr+1 >= len(antennasToPlot)) or |
3270 3395 | ((overlaySpws or overlayBasebands) and spwctr+1 >= len(spwsToPlot))): |
3271 3396 | mytime += 1 |
3272 3397 | if (debug): |
3273 3398 | print("a) xctr=%d, Incrementing mytime to %d" % (xctr, mytime)) |
3274 3399 | if overlayAntennas and xctr+1 == len(antennasToPlot): # second fix for CAS-13568 |
3275 3400 | DrawAntennaNamesForOverlayAntennas(xstartPolLabel, ystartPolLabel, polsToPlot, corr_type, channeldiff, ystartMadLabel, subplotRows, gamp_mad, gamp_std, overlayColors, mysize, ampmarkstyle, markersize, markeredgewidth, msAnt, msFound, antennasToPlot, ampmarkstyle2, xframe, firstFrame, caltableTitle, titlesize, debug=debug) |
3297 3422 | drewAtmosphere = True |
3298 3423 | |
3299 3424 | continue |
3300 3425 | # The following variable allows color legend of UT times to match line plot |
3301 3426 | myUniqueTime = [] |
3302 3427 | if (True): # multiFieldsWithOverlayTime): |
3303 3428 | # support multi-fields with overlay='time' |
3304 3429 | uTPFPS = [] |
3305 3430 | for f in fieldIndicesToPlot: |
3306 3431 | for t in uniqueTimesPerFieldPerSpw[ispwInCalTable][f]: |
3432 + | if tableFormat >= 34: ### added 2024Aug |
3433 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
3434 + | else: ### added 2024Aug |
3435 + | scansToPlotHere = scansToPlot ### added 2024Aug |
3307 3436 | if (sloppyMatch(t, timerangeListTimes, solutionTimeThresholdSeconds, |
3308 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
3309 - | # mytime, scansToPlot, scansForUniqueTimes, # task version |
3437 + | mytime, scansToPlotHere, scansForUniqueTimes, # au version ### modified 2024Aug |
3310 3438 | myprint=debugSloppyMatch |
3311 3439 | )): |
3312 3440 | uTPFPS.append(t) |
3313 3441 | uTPFPS = np.sort(uTPFPS) |
3314 3442 | ctr = 0 |
3315 3443 | for t in uTPFPS: |
3316 3444 | if (debug and False): |
3317 3445 | print("1)checking time %d" % (t)) |
3318 3446 | if (overlayTimes or overlayAntennas): |
3319 3447 | sm = sloppyMatch(uniqueTimes[mytime],times[i],solutionTimeThresholdSeconds,myprint=False) |
3320 3448 | else: |
3449 + | if tableFormat >= 34: ### added 2024Aug |
3450 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
3451 + | else: ### added 2024Aug |
3452 + | scansToPlotHere = scansToPlot ### added 2024Aug |
3321 3453 | sm = sloppyMatch(t, uniqueTimes[mytime], solutionTimeThresholdSeconds, |
3322 - | # mytime, scansToPlot, scansForUniqueTimes, # task version |
3323 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
3454 + | mytime, scansToPlotHere, scansForUniqueTimes, # au version ### modified 2024Aug |
3324 3455 | myprint=debugSloppyMatch |
3325 3456 | ) |
3326 3457 | if (sm): |
3327 3458 | if (debug): |
3328 3459 | print("1)setting myUniqueTime to %d" % (mytime)) |
3329 3460 | myUniqueTime = mytime |
3330 3461 | ctr += 1 |
3331 3462 | if (ctr > len(fieldIndicesToPlot) and bOverlay==False): |
3332 3463 | if (debug): print("multi-field time overlay *************** why are there 2 matches?") |
3333 3464 | # # # # if (ctr == 0): |
3380 3511 | finalTimerangeFlagged = True # 04-Aug-2014 |
3381 3512 | if (debug): |
3382 3513 | print("###### set doneOverlayTime = %s" % (str(doneOverlayTime))) |
3383 3514 | |
3384 3515 | # draw labels |
3385 3516 | # try adding the following 'if' statement on Jun 18, 2013; it works. |
3386 3517 | # if (drewAtmosphere==False or overlayAntennas==False): |
3387 3518 | # Add the 'and not' case to prevent extra atm/fdms shown if one spw's solutions are all flagged |
3388 3519 | if (drewAtmosphere==False or (not overlayAntennas and not allTimesFlaggedOnThisSpw)): |
3389 3520 | if (debug): print("drawOverlayTimeLegends loc 1") |
3390 - | drawOverlayTimeLegends(xframe,firstFrame,xstartTitle,ystartTitle, |
3391 - | caltable,titlesize,fieldIndicesToPlot, |
3392 - | ispwInCalTable,uniqueTimesPerFieldPerSpw, |
3521 + | drawOverlayTimeLegends(xframe, firstFrame, xstartTitle, ystartTitle, |
3522 + | caltable, titlesize, fieldIndicesToPlot, |
3523 + | ispwInCalTable, uniqueTimesPerFieldPerSpw, |
3393 3524 | timerangeListTimes, |
3394 3525 | solutionTimeThresholdSeconds, |
3395 3526 | debugSloppyMatch, |
3396 - | ystartOverlayLegend,debug,mysize, |
3397 - | fieldsToPlot,myUniqueColor, |
3398 - | timeHorizontalSpacing,fieldIndex, |
3527 + | ystartOverlayLegend, debug, mysize, |
3528 + | fieldsToPlot, myUniqueColor, |
3529 + | timeHorizontalSpacing, fieldIndex, |
3399 3530 | overlayColors, |
3400 3531 | antennaVerticalSpacing, overlayAntennas, |
3401 3532 | timerangeList, caltableTitle, mytime, |
3402 - | # scansToPlot, scansForUniqueTimes) # task version |
3403 - | scansToPlotPerSpw[ispw], scansForUniqueTimes) # au version |
3533 + | scansToPlotPerSpw[ispw], scansForUniqueTimes, |
3534 + | uniqueSpwsInCalTable, uniqueTimes) |
3404 3535 | if not LO1 and type(lo1s) == dict: # Fix for SCOPS-4877 |
3405 3536 | LO1 = lo1s[myspw] # Fix for SCOPS-4877 |
3406 3537 | # CAS-8655 |
3407 3538 | newylimits = drawAtmosphereAndFDM(showatm,showtsky,atmString,subplotRows,mysize, |
3408 3539 | TebbSky,TebbSkyImage,plotrange, xaxis,atmchan, |
3409 3540 | atmfreq,transmission,subplotCols,showatmPoints, |
3410 3541 | xframe, channels,LO1,atmchanImage,atmfreqImage, |
3411 3542 | transmissionImage, firstFrame,showfdm,nChannels, |
3412 3543 | tableFormat,originalSpw_casa33, |
3413 3544 | chanFreqGHz_casa33,originalSpw,chanFreqGHz, |
3488 3619 | if (firstTimeMatch == -1): |
3489 3620 | firstTimeMatch = mytime |
3490 3621 | if (debug): |
3491 3622 | print("Setting firstTimeMatch from -1 to %s" % (str(firstTimeMatch))) |
3492 3623 | # The following was needed to support overlay='antenna,time' for showatm for QA2 report (CAS-7820) |
3493 3624 | if (finalTimeMatch == -1 or finalTimeMatch < mytime): |
3494 3625 | if (debug): |
3495 3626 | print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Setting finalTimeMatch from %d to %d" % (finalTimeMatch, mytime)) |
3496 3627 | finalTimeMatch = mytime |
3497 3628 | |
3498 - | ################### Here is the amplitude plotting ############ stopping here Sep 4, 2013 |
3629 + | ################### Here is the amplitude plotting ############ |
3499 3630 | if (yaxis.find('amp')>=0 or yaxis.find('both')>=0 or yaxis.find('ap')>=0) and doneOverlayTime==False: |
3500 3631 | if (overlayBasebands and amplitudeWithPhase): # CAS-6477 |
3501 3632 | if (float(xframe/10) != xframe*0.1 and alreadyPlottedAmp): |
3502 3633 | xframe -= 2 |
3503 3634 | |
3504 3635 | if (debug): |
3505 3636 | print("amp: xctr=%d, xant=%d, myap=%d, mytime=%d(%s), firstTimeMatch=%d, bOverlay=" % (xctr, xant, myap, mytime, utstring(uniqueTimes[mytime],3), firstTimeMatch), bOverlay) |
3506 3637 | if (myap==1): |
3507 3638 | if (overlayTimes == False or mytime==firstTimeMatch): |
3508 3639 | if ((overlaySpws == False and overlayBasebands==False) or spwctr==spwctrFirstToPlot or spwctr>len(spwsToPlot)): |
3706 3837 | markerfacecolor=overlayColors[xctr],markeredgewidth=markeredgewidth) |
3707 3838 | newylimits = recalcYlimits(plotrange,newylimits,gamp[p]) |
3708 3839 | if (overlayAntennas and overlayTimes==False): |
3709 3840 | pb.setp(pdesc, color=overlayColors[xctr]) |
3710 3841 | elif (overlayTimes and overlayAntennas==False): |
3711 3842 | pb.setp(pdesc, color=overlayColors[mytime]) |
3712 3843 | elif (overlayTimes and overlayAntennas): # try to support time,antenna |
3713 3844 | if (debug): |
3714 3845 | print("p=%d, len(fieldsToPlot)=%d, len(timerangeList)=%d" % (p,len(fieldsToPlot),len(timerangeList))) |
3715 3846 | if (len(fieldsToPlot) > 1 or len(timerangeList)>1): |
3716 - | # pb.setp(pdesc, color=overlayColors[myUniqueTime]) |
3717 - | # # # # print("pb.setp: myUniqueTime, overlayColors = ", myUniqueTime, overlayColors[myUniqueTime]) |
3718 3847 | # The third 'or' below is needed if pol='0' is flagged on antenna 0. -- 2012/10/12 |
3719 3848 | if (p==0 or len(polsToPlot)==1 or myUniqueColor==[]): |
3720 3849 | myUniqueColor.append(overlayColors[len(myUniqueColor)]) |
3721 3850 | pb.setp(pdesc, color=myUniqueColor[-1]) |
3722 3851 | else: |
3723 3852 | if (corr_type_string[p] in polsToPlot): |
3724 3853 | # # # # print("pcolor[%d]=%s" % (p,pcolor)) |
3725 3854 | pb.plot(pchannels[p],gamp[p],'%s%s'%(pcolor[p],ampmarkstyle), markersize=markersize,markeredgewidth=markeredgewidth) |
3726 3855 | newylimits = recalcYlimits(plotrange,newylimits,gamp[p]) |
3727 3856 | if (sum(xflag)>0): |
3728 3857 | myxrange = np.max(channels)-np.min(channels) |
3729 3858 | SetNewXLimits([np.min(channels)-myxrange/20, np.max(channels)+myxrange/20],1) |
3730 3859 | # # # # print("amp: Resetting xaxis channel range to counteract flagged data") |
3731 3860 | if (xframe in bottomRowFrames or (xctr+1==len(antennasToPlot) and ispw==spwsToPlot[-1])): |
3732 - | pb.xlabel("Channel", size=mysize) |
3861 + | pb.xlabel("Channels (%d)" % (len(pchannels[p])), size=mysize) ### changed 2024Aug24 |
3733 3862 | elif (xaxis.find('freq')>=0): # amp |
3734 3863 | if (bOverlay): |
3735 3864 | myxrange = np.abs(xfrequencies[0]-xfrequencies[-1]) |
3736 3865 | try: |
3737 3866 | xrange2 = np.abs(xfrequencies2[0]-xfrequencies2[-1]) |
3738 3867 | except: |
3739 3868 | print("No amp data found in second solution. Try increasing the solutionTimeThresholdSeconds above %.0f." % (solutionTimeThresholdSeconds)) |
3740 3869 | print("If this doesn't work, email the developer (%s)." % (developerEmail)) |
3741 3870 | return() |
3742 3871 | |
3760 3889 | pb.plot(pfrequencies2[p], gamp2[p], '%s%s'%(p2color[p],ampmarkstyle), linewidth=width2, markersize=markersize,markeredgewidth=markeredgewidth) |
3761 3890 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp2[p], sideband,plotrange,xchannels2,chanrangePercent=chanrangePercent) |
3762 3891 | for p in range(nPolarizations): |
3763 3892 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
3764 3893 | pb.plot(pfrequencies[p], gamp[p], '%s%s'%(pcolor[p],ampmarkstyle), linewidth=width1, markersize=markersize,markeredgewidth=markeredgewidth) |
3765 3894 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp[p], sideband,plotrange,xchannels,chanrangePercent=chanrangePercent) |
3766 3895 | else: |
3767 3896 | width1 = 1 |
3768 3897 | width2 = 2 # Just enough to distinguish one line from the other. |
3769 3898 | # solutions may be different level of smoothing, so plot highest rms first |
3770 - | if (MAD(gamp[0]) < MAD(gamp2[0])): |
3899 + | if madOfDiff(gamp[0]) < madOfDiff(gamp2[0]): # and firstPlot != 1): # only au version has this parameter |
3900 + | # if (MAD(gamp[0]) < MAD(gamp2[0])): |
3771 3901 | for p in range(nPolarizations): |
3772 3902 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
3773 3903 | pb.plot(pfrequencies2[p], gamp2[p], '%s%s'%(p2color[p],ampmarkstyle), linewidth=width1, markersize=markersize,markeredgewidth=markeredgewidth) |
3774 3904 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp2[p], sideband,plotrange,xchannels2,chanrangePercent=chanrangePercent) |
3775 3905 | for p in range(nPolarizations): |
3776 3906 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
3777 3907 | pb.plot(pfrequencies[p], gamp[p], '%s%s'%(pcolor[p],ampmarkstyle), linewidth=width2, markersize=markersize,markeredgewidth=markeredgewidth) |
3778 3908 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp[p], sideband,plotrange,xchannels,chanrangePercent=chanrangePercent) |
3779 3909 | else: |
3910 + | # plot first solution first |
3780 3911 | for p in range(nPolarizations): |
3781 3912 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
3782 3913 | pb.plot(pfrequencies[p], gamp[p], '%s%s'%(pcolor[p],ampmarkstyle), linewidth=width2, markersize=markersize,markeredgewidth=markeredgewidth) |
3783 3914 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp[p], sideband,plotrange,xchannels,chanrangePercent=chanrangePercent) |
3784 3915 | for p in range(nPolarizations): |
3785 3916 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
3786 3917 | pb.plot(pfrequencies2[p], gamp2[p], '%s%s'%(p2color[p],ampmarkstyle), linewidth=width1, markersize=markersize,markeredgewidth=markeredgewidth) |
3787 3918 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp2[p], sideband,plotrange,xchannels2,chanrangePercent=chanrangePercent) |
3788 3919 | # must set new limits after plotting 'amp' |
3789 3920 | if (zoom=='intersect'): |
3803 3934 | SetLimits(plotrange, chanrange, newylimits, channels, frequencies, |
3804 3935 | pfrequencies, ampMin, ampMax, xaxis, pxl, chanrangeSetXrange, |
3805 3936 | chanrangePercent) |
3806 3937 | else: |
3807 3938 | SetLimits(plotrange, chanrange, newylimits, channels, frequencies2, |
3808 3939 | pfrequencies2, ampMin, ampMax, xaxis, pxl, chanrangeSetXrange, |
3809 3940 | chanrangePercent) |
3810 3941 | # draw polarization and spw labels |
3811 3942 | if (xframe == firstFrame): |
3812 3943 | # draw title including caltable name |
3813 - | caltableList = 'c1 = ' + caltable + ', c2 = ' + caltable2 # + ' (%s)'%(utstring(uniqueTimes2[mytime],3)) |
3944 + | caltableList = 'c1=' + caltable + ', c2=' + caltable2 # + ' (%s)'%(utstring(uniqueTimes2[mytime],3)) |
3814 3945 | pb.text(xstartTitle, ystartTitle, caltableList, size=titlesize, |
3815 3946 | color='k', transform=pb.gcf().transFigure) |
3816 3947 | elif (bpolyOverlay): |
3817 3948 | if (debug): |
3818 3949 | print("in bpolyOverlay **********************************") |
3819 3950 | matches1 = [] |
3951 + | if tableFormat >= 34: ### added 2024Aug |
3952 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
3953 + | else: ### added 2024Aug |
3954 + | scansToPlotHere = scansToPlot ### added 2024Aug |
3820 3955 | for tbp in range(len(timesBP)): |
3821 3956 | if (sloppyMatch(uniqueTimes[mytime], timesBP[tbp], solutionTimeThresholdSeconds, |
3822 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
3823 - | # mytime, scansToPlot, scansForUniqueTimes, # task version |
3957 + | mytime, scansToPlotHere, scansForUniqueTimes, # au version ### modified 2024Aug |
3824 3958 | myprint=debugSloppyMatch |
3825 3959 | )): |
3826 3960 | matches1.append(tbp) |
3827 3961 | matches1 = np.array(matches1) |
3828 3962 | if (len(matches1) < 1): |
3829 3963 | print("No time match found between %.1f and %s" % (uniqueTimes[mytime], str(timesBP))) |
3830 3964 | print("If you are sure the solutions correspond to the same data, you can set solutionTimeThresholdSeconds>=%.0f" % (1+np.ceil(np.abs(timesBP[0]-uniqueTimes[mytime])))) |
3831 3965 | return() |
3832 3966 | matches2 = np.where(xant == np.array(antennasBP))[0] |
3833 3967 | if (len(matches2) < 1): |
3985 4119 | pb.setp(pdesc, color=overlayColors[mytime]) |
3986 4120 | elif (overlayTimes and overlayAntennas): # try to support antenna,time |
3987 4121 | if (myUniqueTime != []): |
3988 4122 | pb.setp(pdesc, color=overlayColors[myUniqueTime]) |
3989 4123 | # The third 'or' below is needed if pol='0' is flagged on antenna 0. -- 2012/10/12 (original spot) |
3990 4124 | if (p==0 or len(polsToPlot)==1 or myUniqueColor==[]): |
3991 4125 | myUniqueColor.append(overlayColors[len(myUniqueColor)]) |
3992 4126 | if (debug): |
3993 4127 | print("myUniqueColor = %s" % (str(myUniqueColor))) |
3994 4128 | pb.setp(pdesc, color=myUniqueColor[-1]) |
3995 - | else: |
4129 + | elif (overlaySpws): # this elif block was missing prior to 2024Aug28 |
4130 + | mycolor = [xcolor,ycolor][p] |
4131 + | linewidth = 1 |
4132 + | pdesc = pb.plot(pfrequencies[p], gamp[p], '%s'%(ampmarkstyles[0]), lw=linewidth, color=mycolor, |
4133 + | markersize=markersize,markeredgewidth=markeredgewidth) |
4134 + | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp[p], sideband, |
4135 + | plotrange,xchannels,debug,-18,chanrangePercent) |
4136 + | else: # show unflagged solutions, no overlay |
3996 4137 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
3997 4138 | # since there is no overlay, don't use dashed line, so zero ------v |
3998 - | pb.plot(pfrequencies[p], gamp[p], '%s%s'%(pcolor[p],ampmarkstyles[0]),markersize=markersize,markeredgewidth=markeredgewidth) |
4139 + | pdesc = pb.plot(pfrequencies[p], gamp[p], '%s%s'%(pcolor[p],ampmarkstyles[0]),markersize=markersize,markeredgewidth=markeredgewidth) |
3999 4140 | newylimits = recalcYlimitsFreq(chanrange, newylimits, gamp[p], sideband,plotrange,xchannels,chanrangePercent=chanrangePercent) |
4000 4141 | # print("newylimits for amp = ", newylimits) |
4001 4142 | # if (debug): print("finished 'for' loop") |
4002 4143 | if (sum(xflag)>0): |
4003 4144 | # print("amp: Resetting xaxis frequency range to counteract flagged data") |
4004 4145 | myxrange = np.max(frequencies)-np.min(frequencies) |
4005 4146 | SetNewXLimits([np.min(frequencies)-0.15*myxrange, np.max(frequencies)+0.15*myxrange],4) |
4006 4147 | |
4007 4148 | if (1==1 or (xframe in bottomRowFrames) or (xctr+1==len(antennasToPlot) and ispw==spwsToPlot[-1])): |
4008 4149 | # use 1==1 because spw might change between top row and bottom row of frames |
4060 4201 | xlim = pb.xlim() |
4061 4202 | ylim = pb.ylim() |
4062 4203 | ResizeFontsSetGrid(adesc,mysize) |
4063 4204 | pb.ylabel(yAmplitudeLabel, size=mysize) |
4064 4205 | pb.subplots_adjust(hspace=myhspace, wspace=mywspace) |
4065 4206 | myxrange = xlim[1]-xlim[0] |
4066 4207 | yrange = ylim[1]-ylim[0] |
4067 4208 | if (debug): print(("amp: ylim, yrange = ", ylim, yrange)) |
4068 4209 | if (overlayAntennas == False and overlayTimes == False and bOverlay == False and |
4069 4210 | ((overlaySpws == False and overlayBasebands == False) or spwctr==spwctrFirstToPlot)): |
4070 - | # draw polarization labels for no overlay |
4211 + | # draw polarization labels for no overlay, or overlaySpws/overlayBasebands |
4071 4212 | x0 = xstartPolLabel |
4072 4213 | y0 = ystartPolLabel |
4073 4214 | for p in range(nPolarizations): |
4074 4215 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
4075 - | pb.text(x0, y0-subplotRows*p*0.03, corrTypeToString(corr_type[p]), |
4076 - | color=pcolor[p],size=mysize, transform=pb.gca().transAxes) |
4077 - | if (channeldiff > 0): |
4078 - | pb.text(x0, ystartMadLabel-0.03*subplotRows*p, |
4079 - | corrTypeToString(corr_type[p])+' MAD = %.4f, St.Dev = %.4f'%(gamp_mad[p]['mad'],gamp_std[p]['std']), |
4080 - | color=pcolor[p],size=mysize, transform=pb.gca().transAxes) |
4216 + | if spwctr==spwctrFirstToPlot or (not overlaySpws and not overlayBasebands): |
4217 + | # no need to plot it more than once in the same position |
4218 + | pb.text(x0, y0-subplotRows*p*0.03, corrTypeToString(corr_type[p]), |
4219 + | color=pcolor[p],size=mysize, transform=pb.gca().transAxes) |
4220 + | if (channeldiff > 0): |
4221 + | pb.text(x0, ystartMadLabel-0.03*subplotRows*p, |
4222 + | corrTypeToString(corr_type[p])+' MAD = %.4f, St.Dev = %.4f'%(gamp_mad[p]['mad'],gamp_std[p]['std']), |
4223 + | color=pcolor[p],size=mysize, transform=pb.gca().transAxes) |
4081 4224 | if (xframe == firstFrame): |
4082 4225 | # draw title including caltable name |
4083 4226 | caltableList = caltableTitle |
4084 4227 | if (bpolyOverlay): |
4085 4228 | caltableList += ', ' + caltable2 + ' (degamp=%d, degphase=%d)'%(nPolyAmp[index]-1,nPolyPhase[index]-1) |
4086 4229 | if (bpolyOverlay2): |
4087 4230 | caltableList += ', ' + caltable3 + ' (degamp=%d, degphase=%d)'%(nPolyAmp2[index]-1,nPolyPhase2[index]-1) |
4088 4231 | pb.text(xstartTitle, ystartTitle, caltableList, size=titlesize, |
4089 4232 | color='k', transform=pb.gcf().transFigure) |
4090 4233 | |
4091 4234 | elif (overlayAntennas==True and xant==antennasToPlot[-1] and bOverlay == False # ): |
4092 - | and overlayTimes==False): # try to support antenna,time avoid antenna labels 'phase' |
4093 - | # We do this last, because by then, the limits will be stable. |
4094 - | if (debug): print("overlayAntennas=True") |
4095 - | x0 = xstartPolLabel |
4096 - | y0 = ystartPolLabel |
4097 - | # draw polarization labels |
4098 - | if (debug): print("1) overlayAntennas=True") |
4099 - | if (corrTypeToString(corr_type[0]) in polsToPlot): |
4100 - | if (channeldiff > 0): |
4101 - | pb.text(x0, ystartMadLabel-0.03*subplotRows*0, |
4102 - | corrTypeToString(corr_type[0])+' MAD = %.4f, St.Dev = %.4f'%(gamp_mad[0]['mad'],gamp_std[0]['std']), |
4103 - | color=overlayColors[0],size=mysize, transform=pb.gca().transAxes) |
4104 - | if (ampmarkstyle.find('-')>=0): |
4105 - | pb.text(x0, y0, corrTypeToString(corr_type[0])+' solid', color=overlayColors[0],size=mysize, |
4106 - | transform=pb.gca().transAxes) |
4107 - | else: |
4108 - | pb.text(x0+0.02, y0, corrTypeToString(corr_type[0]), color=overlayColors[0],size=mysize, |
4109 - | transform=pb.gca().transAxes) |
4110 - | pdesc = pb.plot([x0-0.01], [y0], '%sk'%ampmarkstyle, markersize=markersize, |
4111 - | scalex=False,scaley=False, transform=pb.gca().transAxes,markeredgewidth=markeredgewidth) |
4112 - | if (debug): print("2) overlayAntennas=True") |
4113 - | if (len(corr_type) > 1): |
4114 - | if (corrTypeToString(corr_type[1]) in polsToPlot): |
4115 - | if (channeldiff > 0): |
4116 - | pb.text(x0, ystartMadLabel-0.03*subplotRows*1, |
4117 - | corrTypeToString(corr_type[1])+' MAD = %.4f, St.Dev = %.4f'%(gamp_mad[1]['mad'],gamp_std[1]['std']), |
4118 - | color=overlayColors[0],size=mysize, transform=pb.gca().transAxes) |
4119 - | if (ampmarkstyle2.find('--')>=0): |
4120 - | pb.text(x0, y0-0.03*subplotRows, corrTypeToString(corr_type[1])+' dashed', |
4121 - | color=overlayColors[0],size=mysize, transform=pb.gca().transAxes) |
4122 - | else: |
4123 - | pb.text(x0+0.02, y0-0.03*subplotRows, corrTypeToString(corr_type[1]), |
4124 - | color=overlayColors[0],size=mysize, transform=pb.gca().transAxes) |
4125 - | pdesc = pb.plot([x0-0.01], [y0-0.03*subplotRows], '%sk'%ampmarkstyle2, |
4126 - | markersize=markersize, scalex=False,scaley=False,markeredgewidth=markeredgewidth) |
4127 - | if (debug): print("3) overlayAntennas=True") |
4128 - | if (xframe == firstFrame): |
4129 - | # draw title including caltable name |
4130 - | if (debug): print("4) overlayAntennas=True") |
4131 - | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, color='k', |
4132 - | transform=pb.gcf().transFigure) |
4133 - | if (debug): print("5) overlayAntennas=True") |
4134 - | DrawAntennaNames(msAnt, antennasToPlot, msFound, mysize, overlayColors) |
4135 - | if (debug): print("6) overlayAntennas=True") |
4235 + | and overlayTimes==False): # try to support antenna,time avoid antenna labels 'amp' |
4236 + | # We do this last, because by then, the limits will be stable. |
4237 + | if (debug): print("overlayAntennas=True") |
4238 + | DrawAntennaNamesForOverlayAntennas(xstartPolLabel, ystartPolLabel, polsToPlot, corr_type, channeldiff, ystartMadLabel, subplotRows, gamp_mad, gamp_std, overlayColors, mysize, ampmarkstyle, markersize, markeredgewidth, msAnt, msFound, antennasToPlot, ampmarkstyle2, xframe, firstFrame, caltableTitle, titlesize) |
4136 4239 | elif (overlayTimes==True and bOverlay == False |
4137 4240 | and overlayAntennas==False): # try to support antenna,time |
4138 4241 | doneOverlayTime = True # assumed until proven otherwise in the 'for' loop |
4139 4242 | for f in fieldIndicesToPlot: |
4140 4243 | if (len(uniqueTimesPerFieldPerSpw[ispwInCalTable][f]) > 0): |
4244 + | # this spw/field combination had some entries plotted |
4141 4245 | if ((uniqueTimes[mytime] < uniqueTimesPerFieldPerSpw[ispwInCalTable][f][-1]-solutionTimeThresholdSeconds) and |
4142 - | scansForUniqueTimes[mytime] != scansForUniqueTimes[-1] and # fix for CAS-14096 |
4143 4246 | (uniqueTimes[mytime] < timerangeListTimes[-1])): |
4144 - | if (debug): |
4145 - | print("-----------Not done because %.0f < %.0f-%d for fieldIndex=%d and <%.0f" % (uniqueTimes[mytime], uniqueTimesPerFieldPerSpw[ispwInCalTable][f][-1], solutionTimeThresholdSeconds, f, timerangeListTimes[-1])) |
4146 - | print("-----------ispwInCalTable=%d, mytime=%d, len(uniqueTimes) = %d" % (ispwInCalTable, mytime, len(uniqueTimes))) |
4147 - | doneOverlayTime = False |
4247 + | if tableFormat >= 34: #### added 2024Aug |
4248 + | if scansForUniqueTimes[mytime] not in [scansForUniqueTimes[-1],scansToPlot[-1]] and mytime != timerangeList[-1]: # fix for CAS-14096: if we are on the final time then we are done; and add support for specifying a limited set of timeranges via the timeranges parameter ### added 2024Aug |
4249 + | if (debug): |
4250 + | print("-----------Not done because %.0f < %.0f-%d for fieldIndex=%d and <%.0f and scan%d not in scan[%d,%d] and t%d != t%d" % (uniqueTimes[mytime], uniqueTimesPerFieldPerSpw[ispwInCalTable][f][-1], solutionTimeThresholdSeconds, f, timerangeListTimes[-1],scansForUniqueTimes[mytime],scansForUniqueTimes[-1],scansToPlot[-1],mytime,timerangeList[-1])) |
4251 + | print("-----------ispwInCalTable=%d, mytime=%d, len(uniqueTimes) = %d" % (ispwInCalTable, mytime, len(uniqueTimes))) |
4252 + | doneOverlayTime = False |
4253 + | else: #### added 2024Aug |
4254 + | doneOverlayTime = False #### added 2024Aug |
4148 4255 | if (debug): |
4149 - | print("------doneOverlayTime = %s" % (str(doneOverlayTime))) |
4256 + | print("------doneOverlayTime = %s on mytime %d" % (str(doneOverlayTime),mytime)) |
4150 4257 | if (doneOverlayTime): |
4151 4258 | # either it is the last time of any times in solution, or the last time in the list of times to plot |
4152 4259 | if (debug): |
4153 4260 | print("*** on last time = %d for last fieldIndex %d or %d>=%d" % (mytime,fieldIndex,mytime,timerangeList[-1])) |
4154 4261 | mytime = nUniqueTimes-1 |
4155 4262 | # We do this last, because by then, the limits will be broad enought and stable. |
4156 4263 | # draw polarization labels |
4157 4264 | DrawPolarizationLabelsForOverlayTime(xstartPolLabel,ystartPolLabel,corr_type,polsToPlot, |
4158 4265 | channeldiff,ystartMadLabel,subplotRows,gamp_mad,mysize, |
4159 4266 | ampmarkstyle,markersize,ampmarkstyle2, gamp_std) |
4160 4267 | if (xframe == firstFrame): |
4161 4268 | # draw title including caltable name |
4162 4269 | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, |
4163 4270 | color='k', transform=pb.gcf().transFigure) |
4271 + | if tableFormat >= 34: ### added 2024Aug |
4272 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
4273 + | else: ### added 2024Aug |
4274 + | scansToPlotHere = scansToPlot ### added 2024Aug |
4164 4275 | if (debug): print("drawOverlayTimeLegends loc 2") |
4165 - | drawOverlayTimeLegends(xframe,firstFrame,xstartTitle,ystartTitle, |
4166 - | caltable,titlesize,fieldIndicesToPlot, |
4167 - | ispwInCalTable,uniqueTimesPerFieldPerSpw, |
4276 + | drawOverlayTimeLegends(xframe, firstFrame, xstartTitle, ystartTitle, |
4277 + | caltable, titlesize, fieldIndicesToPlot, |
4278 + | ispwInCalTable, uniqueTimesPerFieldPerSpw, |
4168 4279 | timerangeListTimes, solutionTimeThresholdSeconds, |
4169 - | debugSloppyMatch,ystartOverlayLegend,debug,mysize, |
4170 - | fieldsToPlot,myUniqueColor,timeHorizontalSpacing, |
4280 + | debugSloppyMatch, ystartOverlayLegend, debug, mysize, |
4281 + | fieldsToPlot, myUniqueColor,timeHorizontalSpacing, |
4171 4282 | fieldIndex, overlayColors, antennaVerticalSpacing, |
4172 4283 | overlayAntennas, timerangeList, caltableTitle, |
4173 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes) |
4174 - | # mytime, scansToPlot, scansForUniqueTimes) # task version |
4175 - | if (debug): print("done drawOverlayTimeLegends loc 2") |
4284 + | mytime, scansToPlotHere, scansForUniqueTimes, |
4285 + | uniqueSpwsInCalTable, uniqueTimes) ### modified 2024Aug |
4176 4286 | elif (overlayAntennas and overlayTimes): # Oct 23, 2012 |
4177 4287 | # This will only happen for overlay='antenna,time' |
4178 4288 | if (xframe == firstFrame and mytime == firstTimeMatch and xctr==firstUnflaggedAntennaToPlot and bOverlay==False): # bug fix on 2015-08-19 for CAS-7820 |
4179 4289 | # draw title including caltable name |
4180 4290 | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, color='k', |
4181 4291 | transform=pb.gcf().transFigure) |
4182 4292 | DrawBottomLegendPageCoords(msName, uniqueTimes[mytime], mysize, figfile) |
4183 4293 | # Adding the following 'for' loop on Mar 13, 2013 to support the case of |
4184 4294 | # single time range with overlay='antenna,time' |
4185 4295 | if (xant==antennasToPlot[-1]): |
4186 4296 | doneOverlayTime = True # assumed until proven otherwise in the 'for' loop |
4187 4297 | for f in fieldIndicesToPlot: |
4188 4298 | if (len(uniqueTimesPerFieldPerSpw[ispwInCalTable][f]) > 0): |
4189 4299 | if ((uniqueTimes[mytime] < uniqueTimesPerFieldPerSpw[ispwInCalTable][f][-1]-solutionTimeThresholdSeconds) and |
4190 4300 | (uniqueTimes[mytime] < timerangeListTimes[-1])): |
4191 - | if (debug): |
4301 + | if tableFormat >= 34: #### added 2024Aug |
4302 + | # if (scansForUniqueTimes[mytime] != scansForUniqueTimes[-1]) #### added 2024Aug (fix for CAS-14096): if we are on the final time then we are done |
4303 + | if scansForUniqueTimes[mytime] != scansForUniqueTimes[-1] and mytime != timerangeList[-1]: # fix for CAS-14096: if we are on the final time then we are done; and add support for specifying a limited set of timeranges via the timeranges parameter ### added 2024Aug |
4304 + | if (debug): |
4192 4305 | print("-----------Not done because %.0f < %.0f-%d for fieldIndex=%d and <%.0f" % (uniqueTimes[mytime], uniqueTimesPerFieldPerSpw[ispwInCalTable][f][-1], solutionTimeThresholdSeconds, f, timerangeListTimes[-1])) |
4193 4306 | print("-----------ispwInCalTable=%d, mytime=%d, len(uniqueTimes) = %d" % (ispwInCalTable, mytime, len(uniqueTimes))) |
4194 - | doneOverlayTime = False |
4307 + | doneOverlayTime = False |
4308 + | else: #### added 2024Aug |
4309 + | doneOverlayTime = False #### added 2024Aug |
4195 4310 | if (doneOverlayTime): |
4196 4311 | # This is necessary for the case that no antennas were flagged for the single timerange selected |
4197 4312 | if (debug): print("drawOverlayTimeLegends loc 3") |
4198 4313 | drawOverlayTimeLegends(xframe,firstFrame,xstartTitle,ystartTitle,caltable,titlesize, |
4199 4314 | fieldIndicesToPlot,ispwInCalTable,uniqueTimesPerFieldPerSpw, |
4200 4315 | timerangeListTimes, solutionTimeThresholdSeconds, |
4201 4316 | debugSloppyMatch,ystartOverlayLegend,debug,mysize, |
4202 4317 | fieldsToPlot,myUniqueColor,timeHorizontalSpacing, |
4203 4318 | fieldIndex, overlayColors, antennaVerticalSpacing, |
4204 4319 | overlayAntennas, timerangeList, caltableTitle, |
4205 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes) |
4206 - | # mytime, scansToPlot, scansForUniqueTimes) # task version |
4320 + | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, |
4321 + | uniqueSpwsInCalTable, uniqueTimes) |
4207 4322 | |
4323 + | # endif / elif / elif / elif |
4208 4324 | |
4209 4325 | if (debug): print("####### 2nd place") |
4210 4326 | # Here is 2nd place where we eliminate any white space on the right and left edge of the plots: 'amp' |
4211 4327 | # |
4212 4328 | if (abs(plotrange[2]) > 0 or abs(plotrange[3]) > 0): |
4213 4329 | SetNewYLimits([plotrange[2],plotrange[3]]) |
4214 4330 | if (plotrange[0]==0 and plotrange[1]==0): |
4215 4331 | if (xaxis.find('chan')>=0): |
4216 4332 | SetNewXLimits([channels[0],channels[-1]],9) |
4217 4333 | else: |
4234 4350 | # I need the following line for chanrange to work |
4235 4351 | if (chanrange[0] != 0 or chanrange[1] != 0 or chanrangePercent != None): |
4236 4352 | SetLimits(plotrange, chanrange, newylimits, channels, frequencies, pfrequencies, |
4237 4353 | ampMin, ampMax, xaxis,pxl, chanrangeSetXrange, |
4238 4354 | chanrangePercent) |
4239 4355 | |
4240 4356 | # Finally, draw the atmosphere and FDM windows, if requested. 'amp' |
4241 4357 | if ((overlayAntennas==False and overlayTimes==False) or |
4242 4358 | (overlayAntennas==True and overlayTimes==False and xant==antennasToPlot[-1]) or |
4243 4359 | (overlayTimes==True and overlayAntennas==False and doneOverlayTime) or |
4244 - | # # (xant==antennasToPlot[-1] and doneOverlayTime) # support showatm with overlay='antenna,time' |
4245 4360 | (overlayTimes and overlayAntennas and # Aug 5, 2013 |
4246 4361 | xant==antennasToPlot[-1] and doneOverlayTime and mytime==finalTimeMatch # 2015-08-19 for CAS-7820 |
4247 4362 | and not drewAtmosphere) # added on 2014-12-04 to support case of a flagged antenna CAS-7187 |
4248 4363 | ): |
4249 4364 | if ((showatm or showtsky) and len(atmString) > 0): |
4250 4365 | DrawAtmosphere(showatm, showtsky, subplotRows, atmString, |
4251 4366 | mysize, TebbSky, plotrange, xaxis, atmchan, |
4252 4367 | atmfreq, transmission, subplotCols, |
4253 4368 | showatmPoints=showatmPoints, xframe=xframe, |
4254 4369 | channels=channels,mylineno=lineNumber(), |
4295 4410 | for p in range(nPolarizations): |
4296 4411 | if (corrTypeToString(corr_type[0]) in polsToPlot): |
4297 4412 | pb.text(x0+0.1, y0-p*0.03*subplotRows, corrTypeToString(corr_type[p]), color=p2color[p], |
4298 4413 | size=mysize,transform=pb.gca().transAxes) |
4299 4414 | if (bpolyOverlay2): |
4300 4415 | for p in range(nPolarizations): |
4301 4416 | if (corrTypeToString(corr_type[0]) in polsToPlot): |
4302 4417 | pb.text(x0+0.2, y0-p*0.03*subplotRows, corrTypeToString(corr_type[p]), |
4303 4418 | color=p3color[p], size=mysize,transform=pb.gca().transAxes) |
4304 4419 | |
4305 - | # # # # if (xframe == 111 and amplitudeWithPhase): |
4306 4420 | myIndexTime = uniqueTimesPerFieldPerSpw[ispwInCalTable][fieldIndex][-1] |
4307 4421 | if (debug): print("running sloppyMatch") |
4422 + | if tableFormat >= 34: ### added 2024Aug |
4423 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
4424 + | else: ### added 2024Aug |
4425 + | scansToPlotHere = scansToPlot ### added 2024Aug |
4308 4426 | matched,mymatch = sloppyMatch(myIndexTime,uniqueTimes,solutionTimeThresholdSeconds, |
4309 - | mytime, scansToPlotPerSpw[ispw], # add PerSpw[ispw] on 2014-04-05 |
4427 + | mytime, scansToPlotHere, ### modified 2024Aug |
4310 4428 | scansForUniqueTimes, |
4311 4429 | whichone=True,myprint=debug) |
4312 4430 | if (debug): |
4313 4431 | print("1)done sloppyMatch, mytime=%d, scansForUniqueTimes=%s" % (mytime,str(scansForUniqueTimes))) |
4314 4432 | print("ispw=%d" % (ispw)) |
4315 4433 | print("len(scansToPlotPerSpw)=%d" % (len(scansToPlotPerSpw))) |
4316 4434 | if (matched == False and scansForUniqueTimes[mytime] in scansToPlotPerSpw[ispw]): |
4317 4435 | print("---------- 1) Did not find %f in %s" % (myIndexTime,str(uniqueTimes))) |
4318 4436 | print("Try re-running with a smaller solutionTimeThresholdSeconds (currently %f)" % (solutionTimeThresholdSeconds)) |
4319 4437 | return |
4320 4438 | else: |
4321 4439 | # we are on the final time to be plotted |
4322 4440 | if (debug): print("on the final time") |
4323 4441 | mytimeTest = mytime==nUniqueTimes-1 # mytime==myIndexTime # mytime==mymatch |
4324 4442 | if ((xframe == 111 and amplitudeWithPhase) or |
4325 4443 | # Following case is needed to make subplot=11 to work for: try to support overlay='antenna,time' |
4326 4444 | (xframe == lastFrame and overlayTimes and overlayAntennas and |
4327 4445 | xctr+1==len(antennasToPlot) and |
4328 - | # # mytime+1==len(uniqueTimes) and # this worked for nspw <= 4 |
4329 4446 | mytimeTest and |
4330 4447 | spwctr<len(spwsToPlot))): # removed +1 from spwctr+1 on 2014-04-05 to match au |
4331 4448 | if (debug): |
4332 4449 | print("xframe=%d == lastFrame=%d, amplitudeWithPhase=%s" % (xframe, lastFrame, str(amplitudeWithPhase))) |
4333 4450 | print("xctr+1=%d == len(antennasToPlot)=%d" % (xctr+1,len(antennasToPlot))) |
4334 4451 | print("mytime+1=%d == len(uniqueTimes)=%d" % (mytime+1,len(uniqueTimes))) |
4335 4452 | print("spwctr+1=%d < len(spwsToPlot)=%d" % (spwctr+1,len(spwsToPlot))) |
4336 4453 | if (len(figfile) > 0): |
4337 4454 | plotfiles.append(makeplot(figfile,msFound,msAnt, |
4338 4455 | overlayAntennas,pages,pagectr, |
4339 4456 | density,interactive,antennasToPlot, |
4340 4457 | spwsToPlot,overlayTimes,overlayBasebands, |
4341 4458 | 3,xant,ispw,subplot,resample, |
4342 4459 | debug,figfileSequential,figfileNumber)) |
4343 4460 | figfileNumber += 1 |
4344 4461 | |
4345 4462 | donetime = time.time() |
4346 4463 | drewAtmosphere = False # needed for CAS-7187 (subplot=11) |
4347 4464 | if (interactive): |
4348 4465 | pb.draw() |
4349 - | # # # # myinput = raw_input(":(%.1f sec) Press return for next page (b for backwards, q to quit): "%(donetime-mytimestamp)) |
4350 4466 | myinput = input("Press return for next page (b for backwards, q to quit): ") |
4351 4467 | else: |
4352 4468 | myinput = '' |
4353 4469 | skippingSpwMessageSent = 0 |
4354 4470 | mytimestamp = time.time() |
4355 4471 | if (myinput.find('q') >= 0): |
4356 4472 | showFinalMessage(overlayAntennas, solutionTimeSpread, nUniqueTimes) |
4357 4473 | return() |
4358 4474 | if (myinput.find('b') >= 0): |
4359 4475 | if (pagectr > 0): |
4416 4532 | if (bOverlay): |
4417 4533 | pchannels2 = [xchannels2,ychannels2] # this is necessary because np.diff reduces nchan by 1 |
4418 4534 | pfrequencies2 = [xfrequencies2,yfrequencies2] # this is necessary because np.diff reduces nchan by 1 |
4419 4535 | if (overlayTimes == False or mytime==firstTimeMatch): |
4420 4536 | if ((overlaySpws == False and overlayBasebands==False) or spwctr==spwctrFirstToPlot or |
4421 4537 | spwctr>spwsToPlot[-1] or |
4422 4538 | (overlayBasebands and amplitudeWithPhase)): # CAS-6477 |
4423 4539 | if (overlayAntennas==False or xctr==firstUnflaggedAntennaToPlot |
4424 4540 | or xctr>antennasToPlot[-1]): # 2012-05-24, to fix the case where all ants flagged on one timerange |
4425 4541 | xframe += 1 |
4426 - | # # # # print("u) incrementing xframe to %d" % xframe) |
4542 + | if debug: |
4543 + | print("u) incrementing xframe to %d" % xframe) |
4427 4544 | myUniqueColor = [] |
4428 4545 | newylimits = [LARGE_POSITIVE, LARGE_NEGATIVE] |
4429 4546 | if (phase != ''): |
4430 4547 | if ((phase[0] != 0 or phase[1] != 0) and amplitudeWithPhase): |
4431 4548 | newylimits = phase |
4432 4549 | if (debug): |
4433 4550 | print("$$$$$$$$$$$$$$$$$$$$$$$ ready to plot phase on xframe %d" % (xframe)) |
4434 4551 | |
4435 4552 | |
4436 4553 | if (previousSubplot != xframe): |
4552 4669 | newylimits = [-minPhaseRange,minPhaseRange] |
4553 4670 | if (phase != ''): |
4554 4671 | if ((phase[0] != 0 or phase[1] != 0) and amplitudeWithPhase): |
4555 4672 | newylimits = phase |
4556 4673 | if (sum(xflag)>0): |
4557 4674 | # # # # print("phase: Resetting xaxis channel range to counteract flagged data") |
4558 4675 | myxrange = np.max(channels)-np.min(channels) |
4559 4676 | SetNewXLimits([np.min(channels)-myxrange/20, |
4560 4677 | np.max(channels)+myxrange/20],14) |
4561 4678 | if (xframe in bottomRowFrames or (xctr+1==len(antennasToPlot) and ispw==spwsToPlot[-1])): |
4562 - | pb.xlabel("Channel", size=mysize) |
4679 + | pb.xlabel("Channels (%d)"%(len(pchannels[p])), size=mysize) ### changed 2024Aug24 |
4563 4680 | elif (xaxis.find('freq')>=0): # 'phase' |
4564 4681 | if (bOverlay): |
4565 4682 | if (debug): |
4683 + | p = 0 ### added 2024Aug to prevent crash due to undefined variable |
4566 4684 | print("Preparing to plot phase from %f-%f for pols: %s" % (xfrequencies[0],xfrequencies[-1],str(polsToPlot))) |
4567 4685 | print("Preparing to plot phase from %f-%f for pols: %s" % (pfrequencies[p][0],pfrequencies[p][-1],str(polsToPlot))) |
4568 4686 | print("Preparing to plot phase from %f-%f for pols: %s" % (pfrequencies2[p][0],pfrequencies2[p][-1],str(polsToPlot))) |
4569 4687 | myxrange = np.abs(xfrequencies[0]-xfrequencies[-1]) |
4570 4688 | try: |
4571 4689 | xrange2 = np.abs(xfrequencies2[0]-xfrequencies2[-1]) |
4572 4690 | except: |
4573 4691 | print("No phase data found in second solution. Try increasing the solutionTimeThresholdSeconds above %.0f." % (solutionTimeThresholdSeconds)) |
4574 4692 | print("If this doesn't work, email the developer (%s)." % (developerEmail)) |
4575 4693 | return() |
4647 4765 | SetLimits(plotrange, chanrange, newylimits, channels, frequencies, |
4648 4766 | pfrequencies, ampMin, ampMax, xaxis,pxl, chanrangeSetXrange, |
4649 4767 | chanrangePercent) |
4650 4768 | else: |
4651 4769 | SetLimits(plotrange, chanrange, newylimits, channels, frequencies2, |
4652 4770 | pfrequencies2, ampMin, ampMax, xaxis,pxl, chanrangeSetXrange, |
4653 4771 | chanrangePercent) |
4654 4772 | # draw polarization and spw labels |
4655 4773 | if (xframe == firstFrame): |
4656 4774 | # draw title including caltable name |
4657 - | caltableList = 'c1 = ' + caltable + ', c2 = ' + caltable2 # + ' (%s)'%(utstring(uniqueTimes2[mytime],3)) |
4775 + | caltableList = 'c1=' + caltable + ', c2=' + caltable2 # + ' (%s)'%(utstring(uniqueTimes2[mytime],3)) |
4658 4776 | pb.text(xstartTitle, ystartTitle, caltableList, size=titlesize, |
4659 4777 | color='k', transform=pb.gcf().transFigure) |
4660 4778 | elif (bpolyOverlay): |
4661 4779 | matches1 = [] |
4780 + | if tableFormat >= 34: ### added 2024Aug |
4781 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
4782 + | else: ### added 2024Aug |
4783 + | scansToPlotHere = scansToPlot ### added 2024Aug |
4662 4784 | for tbp in range(len(timesBP)): |
4663 4785 | if (sloppyMatch(uniqueTimes[mytime], timesBP[tbp], solutionTimeThresholdSeconds, |
4664 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
4665 - | # mytime, scansToPlot, scansForUniqueTimes, # task version |
4786 + | mytime, scansToPlotHere, scansForUniqueTimes, ### modified 2024Aug |
4666 4787 | myprint=debugSloppyMatch)): |
4667 4788 | matches1.append(tbp) |
4668 4789 | matches1 = np.array(matches1) |
4669 4790 | # # # # print("time matches: matches1 = ", matches1) |
4670 4791 | if (len(matches1) < 1): |
4671 4792 | print("No time match found") |
4672 4793 | print("If you are sure the solutions correspond to the same data, you can set solutionTimeThresholdSeconds=%.0f" % (1+np.ceil(np.abs(timesBP[0]-uniqueTimes[mytime])))) |
4673 4794 | return() |
4674 - | # # # # matches1 = np.where(np.floor(uniqueTimes[mytime]) == np.floor(np.array(timesBP)))[0] |
4675 4795 | matches2 = np.where(xant == np.array(antennasBP))[0] |
4676 4796 | if (len(matches2) < 1): |
4677 4797 | print("No antenna match found between %s and %s" % (str(xant), str(antennasBP))) |
4678 4798 | # # # # print("antenna matches: matches2 = ", matches2) |
4679 4799 | |
4680 4800 | if (tableFormat == 33): |
4681 4801 | matches3 = np.where(ispw == np.array(cal_desc_idBP))[0] |
4682 4802 | if (len(matches3) < 1): |
4683 4803 | print("No spw match found: %d not in %s" % (ispw, str(cal_desc_idBP))) |
4684 4804 | else: |
4898 5018 | caltableList = caltableTitle |
4899 5019 | if (bpolyOverlay): |
4900 5020 | caltableList += ', ' + caltable2 + ' (degamp=%d, degphase=%d)'%(nPolyAmp[index]-1,nPolyPhase[index]-1) |
4901 5021 | if (bpolyOverlay2): |
4902 5022 | caltableList += ', ' + caltable3 + ' (degamp=%d, degphase=%d)'%(nPolyAmp2[index]-1,nPolyPhase2[index]-1) |
4903 5023 | pb.text(xstartTitle, ystartTitle, caltableList, size=titlesize, |
4904 5024 | color='k', transform=pb.gcf().transFigure) |
4905 5025 | elif (overlayAntennas==True and xant==antennasToPlot[-1] and bOverlay==False # ): |
4906 5026 | and overlayTimes==False): # try to support antenna,time avoid antenna labels 'phase' |
4907 5027 | # We do this last, because by then, the limits will be stable. |
5028 + | # draw polarization labels for overlayAntennas |
4908 5029 | x0 = xstartPolLabel |
4909 5030 | y0 = ystartPolLabel |
4910 5031 | if (corrTypeToString(corr_type[0]) in polsToPlot): |
4911 5032 | if (channeldiff > 0): |
4912 5033 | pb.text(x0, ystartMadLabel-0.03*subplotRows*p, |
4913 5034 | corrTypeToString(corr_type[p])+' MAD = %.4f, St.Dev = %.4f'%(gphs_mad[p]['mad'],gphs_std[p]['std']), |
4914 5035 | color=overlayColors[0], size=mysize, transform=pb.gca().transAxes) |
4915 5036 | if (phasemarkstyle.find('-')>=0): |
4916 5037 | pb.text(x0, y0-0.03*subplotRows*0, corrTypeToString(corr_type[0])+' solid', color=overlayColors[0], |
4917 5038 | fontsize=mysize, transform=pb.gca().transAxes) |
4937 5058 | if (xframe == firstFrame): |
4938 5059 | # draw title including caltable name |
4939 5060 | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, color='k', |
4940 5061 | transform=pb.gcf().transFigure) |
4941 5062 | DrawAntennaNames(msAnt, antennasToPlot, msFound, mysize, overlayColors) |
4942 5063 | elif (overlayTimes==True and bOverlay == False |
4943 5064 | and overlayAntennas==False): # try to support antenna,time |
4944 5065 | doneOverlayTime = True # assumed until proven otherwise in the 'for' loop |
4945 5066 | for f in fieldIndicesToPlot: |
4946 5067 | if (uniqueTimes[mytime] < uniqueTimesPerFieldPerSpw[ispwInCalTable][f][-1]-solutionTimeThresholdSeconds and |
4947 - | scansForUniqueTimes[mytime] != scansForUniqueTimes[-1] and # fix for CAS-14096 |
4948 5068 | uniqueTimes[mytime] < timerangeListTimes[-1]): |
4949 - | doneOverlayTime = False |
5069 + | if tableFormat >= 34: # added 2024Aug |
5070 + | if scansForUniqueTimes[mytime] != scansForUniqueTimes[-1]: # fix for CAS-14096 |
5071 + | doneOverlayTime = False |
5072 + | else: # added 2024Aug |
5073 + | doneOverlayTime = False # added 2024Aug |
4950 5074 | if (doneOverlayTime): |
4951 5075 | # either it is the last time of any times in solution, or the last time in the list of times to plot |
4952 5076 | mytime = nUniqueTimes-1 |
5077 + | # draw polarization labels for overlayTimes |
4953 5078 | # We do this last, because by then, the limits will be broad enough and stable. |
4954 5079 | x0 = xstartPolLabel |
4955 5080 | y0 = ystartPolLabel |
4956 5081 | if (corrTypeToString(corr_type[0]) in polsToPlot): |
4957 5082 | if (channeldiff > 0): |
4958 5083 | p = 0 |
4959 5084 | pb.text(x0, ystartMadLabel-0.03*subplotRows*p, |
4960 5085 | corrTypeToString(corr_type[p])+' MAD = %.4f'%(gphs_mad[p]['mad']), |
4961 5086 | color='k', size=mysize, transform=pb.gca().transAxes) |
4962 5087 | if (phasemarkstyle.find('-')>=0): |
4988 5113 | transform=pb.gcf().transFigure) |
4989 5114 | if (debug): print("drawOverlayTimeLegends loc 4") |
4990 5115 | drawOverlayTimeLegends(xframe,firstFrame,xstartTitle,ystartTitle, |
4991 5116 | caltable,titlesize,fieldIndicesToPlot, |
4992 5117 | ispwInCalTable,uniqueTimesPerFieldPerSpw, |
4993 5118 | timerangeListTimes, solutionTimeThresholdSeconds, |
4994 5119 | debugSloppyMatch,ystartOverlayLegend,debug,mysize, |
4995 5120 | fieldsToPlot,myUniqueColor,timeHorizontalSpacing, |
4996 5121 | fieldIndex, overlayColors, antennaVerticalSpacing, |
4997 5122 | overlayAntennas, timerangeList, caltableTitle, |
4998 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes) |
4999 - | # mytime, scansToPlot, scansForUniqueTimes) # task version |
5123 + | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, |
5124 + | uniqueSpwsInCalTable, uniqueTimes) |
5000 5125 | |
5001 5126 | elif (overlayAntennas and overlayTimes): # Oct 23, 2012 |
5002 5127 | # This will only happen for: try to support overlay='antenna,time' |
5003 5128 | if (xframe == firstFrame and mytime==0 and xctr==firstUnflaggedAntennaToPlot and bOverlay==False): |
5004 5129 | # draw title including caltable name |
5005 5130 | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, color='k', |
5006 5131 | transform=pb.gcf().transFigure) |
5007 5132 | DrawBottomLegendPageCoords(msName, uniqueTimes[mytime], mysize, figfile) |
5008 5133 | |
5009 5134 | #endif (overlayAntennas == False and overlayTimes == False and bOverlay == False) |
5067 5192 | showtsys=showtsys, Trx=Trx) |
5068 5193 | drewAtmosphere = True |
5069 5194 | |
5070 5195 | if (xaxis.find('freq')>=0 and showfdm and nChannels <= 256): |
5071 5196 | if (tableFormat == 33): |
5072 5197 | showFDM(originalSpw_casa33, chanFreqGHz_casa33, baseband, showBasebandNumber, basebandDict, overlayColors) |
5073 5198 | else: |
5074 5199 | showFDM(originalSpw, chanFreqGHz, baseband, showBasebandNumber, basebandDict, overlayColors) |
5075 5200 | |
5076 5201 | if (bOverlay): |
5077 - | # draw polarization labels |
5202 + | # draw polarization labels for bOverlay |
5078 5203 | x0 = xstartPolLabel |
5079 5204 | y0 = ystartPolLabel |
5080 5205 | for p in range(nPolarizations): |
5081 5206 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
5082 5207 | pb.text(x0, y0-p*0.03*subplotRows, corrTypeToString(corr_type[p])+'-c1', |
5083 5208 | color=pcolor[p],size=mysize,transform=pb.gca().transAxes) |
5084 5209 | pb.text(x0, y0-(p*0.03+0.06)*subplotRows, corrTypeToString(corr_type[p])+'-c2', |
5085 5210 | color=p2color[p],size=mysize, transform=pb.gca().transAxes) |
5086 5211 | if (bpolyOverlay and xaxis.find('freq')>=0): |
5212 + | # draw polarization labels for bpolyOverlay |
5087 5213 | x0 = xstartPolLabel |
5088 5214 | y0 = ystartPolLabel |
5089 5215 | if (xcolor != x2color): |
5090 5216 | for p in range(nPolarizations): |
5091 5217 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
5092 5218 | pb.text(x0+0.1, y0-p*0.03*subplotRows, corrTypeToString(corr_type[p]), color=p2color[p], |
5093 5219 | size=mysize, transform=pb.gca().transAxes) |
5094 5220 | if (bpolyOverlay2): |
5095 5221 | for p in range(nPolarizations): |
5096 5222 | if (corrTypeToString(corr_type[p]) in polsToPlot): |
5102 5228 | redisplay = False |
5103 5229 | |
5104 5230 | if (xframe == lastFrame): |
5105 5231 | if (debug): |
5106 5232 | print("*** mytime+1=%d, nUniqueTimes=%d, timerangeList[-1]=%d, doneOverlayTime=%s" % (mytime+1, nUniqueTimes,timerangeList[-1],doneOverlayTime)) |
5107 5233 | print("*** xant=%d, antennasToPlot[-1]=%d, overlayAntennas=%s, overlayTimes=%s" % (xant,antennasToPlot[-1],overlayAntennas,overlayTimes)) |
5108 5234 | print("*** xframe=%d, lastFrame=%d, xctr=%d, spwctr=%d, len(antennasToPlot)=%d, len(spwsToPlot)=%d" % (xframe,lastFrame,xctr,spwctr,len(antennasToPlot), len(spwsToPlot))) |
5109 5235 | myIndexTime = uniqueTimesPerFieldPerSpw[ispwInCalTable][fieldIndex][-1] |
5110 5236 | if (debug): |
5111 5237 | print("myIndexTime = ", myIndexTime) |
5238 + | if tableFormat >= 34: ### added 2024Aug |
5239 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
5240 + | else: ### added 2024Aug |
5241 + | scansToPlotHere = scansToPlot ### added 2024Aug |
5112 5242 | matched,mymatch = sloppyMatch(myIndexTime,uniqueTimes,solutionTimeThresholdSeconds, |
5113 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
5114 - | # mytime, scansToPlot, scansForUniqueTimes, # task version |
5243 + | mytime, scansToPlotHere, scansForUniqueTimes, |
5115 5244 | whichone=True, myprint=False) |
5116 5245 | if (matched==False and scansForUniqueTimes[mytime] in scansToPlotPerSpw[ispw]): |
5117 5246 | print("---------- 2) Did not find %f within %.1f seconds of anything in %s" % (myIndexTime,solutionTimeThresholdSeconds,str(uniqueTimes))) |
5118 5247 | print("Try re-running with a smaller solutionTimeThresholdSeconds (currently %f)" % (solutionTimeThresholdSeconds)) |
5119 5248 | return |
5120 5249 | else: |
5121 5250 | # we are on the final time to be plotted |
5122 5251 | if (debug): print("on the final time") |
5123 5252 | mytimeTest = mytime==nUniqueTimes-1 |
5124 5253 | if (debug): |
5127 5256 | # old 3.3 cal tables will land here |
5128 5257 | scanTest = False |
5129 5258 | scanTest2 = False |
5130 5259 | else: |
5131 5260 | if (debug): |
5132 5261 | print("ispw=%d len(scansToPlotPerSpw[ispw])=%d mytime=%d, len(scansForUniqueTimes)=%d" % (ispw,len(scansToPlotPerSpw[ispw]),mytime,len(scansForUniqueTimes))) |
5133 5262 | print("scansToPlotPerSpw = ", scansToPlotPerSpw) |
5134 5263 | if (len(scansToPlotPerSpw[ispw]) == 0): |
5135 5264 | scanTest = False |
5136 5265 | else: |
5137 - | scanTest = scansToPlotPerSpw[ispw][-1]==scansForUniqueTimes[mytime] |
5266 + | scanTest = scansToPlotPerSpw[ispw][-1]==scansForUniqueTimes[mytime] and (VisCal != 'SDSKY_PS' or mytime == timerangeList[-1]) ### added second expression 2024Aug27 |
5138 5267 | highestSpwIndexInSpwsToPlotThatHasCurrentScan = \ |
5139 5268 | computeHighestSpwIndexInSpwsToPlotThatHasCurrentScan(spwsToPlot, scansToPlotPerSpw, scansForUniqueTimes[mytime]) |
5140 5269 | if (highestSpwIndexInSpwsToPlotThatHasCurrentScan == -1): |
5141 5270 | scanTest2 = False |
5142 5271 | else: |
5143 5272 | scanTest2 = (spwctr == highestSpwIndexInSpwsToPlotThatHasCurrentScan) |
5144 5273 | if ((overlayAntennas==False and overlayTimes==False and overlaySpws==False and overlayBasebands==False) |
5145 5274 | # either it is the last time of any, or the last time in the list of times to plot |
5146 5275 | or (overlayAntennas==False and overlaySpws==False and overlayBasebands==False and (mytime+1==nUniqueTimes or mytime == timerangeList[-1])) # or mytimeTest)) # removed on July 25,2013 |
5147 5276 | or (xant==antennasToPlot[-1] and overlayAntennas==True and overlayTimes==False and overlaySpws==False and overlayBasebands==False) |
5216 5345 | xant = antennasToPlot[xctr] |
5217 5346 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
5218 5347 | ispw = spwsToPlot[spwctr] |
5219 5348 | # # # # print("Returning to [%d,%d,%d,%d]" % (xctr,spwctr,mytime,myap)) |
5220 5349 | else: |
5221 5350 | pagectr += 1 |
5222 5351 | if (pagectr >= len(pages)): |
5223 5352 | newpage = 1 |
5224 5353 | else: |
5225 5354 | newpage = 0 |
5355 + | if tableFormat >= 34: ### added 2024Aug |
5356 + | scansToPlotHere = scansToPlotPerSpw[ispw] ### added 2024Aug |
5357 + | else: ### added 2024Aug |
5358 + | scansToPlotHere = scansToPlot ### added 2024Aug |
5226 5359 | if (overlayTimes==True and |
5227 5360 | sloppyMatch(uniqueTimesPerFieldPerSpw[ispwInCalTable][fieldIndex][-1], |
5228 5361 | uniqueTimes[mytime],solutionTimeThresholdSeconds, |
5229 - | mytime, scansToPlotPerSpw[ispw], scansForUniqueTimes, # au version |
5230 - | # mytime, scansToPlot, scansForUniqueTimes, # task version |
5362 + | mytime, scansToPlotHere, scansForUniqueTimes, # au version ### modified 2024Aug |
5231 5363 | myprint=debugSloppyMatch)): |
5232 5364 | # be sure to avoid any more loops through mytime which will cause 'b' button to fail |
5233 - | mytime = nUniqueTimes |
5365 + | if VisCal != 'SDSKY_PS': |
5366 + | mytime = nUniqueTimes |
5234 5367 | else: |
5235 5368 | if (debug): |
5236 5369 | print(">>>>>>>>>>> Not going to new page, uniqueTimes[mytime]=%.8f, uniqueTimesPerFieldPerSpw[ispwInCalTable=%d][fieldIndex=%d][-1]=%.8f" % (uniqueTimes[mytime], ispwInCalTable, fieldIndex, uniqueTimesPerFieldPerSpw[ispwInCalTable][fieldIndex][-1])) |
5237 5370 | print("spwctr=%d ?== (len(spwsToPlot)-1)=%d, spwsToPlot=" % (spwctr,len(spwsToPlot)-1),spwsToPlot) |
5238 5371 | print("test1: %s" % (overlayAntennas==False and overlayTimes==False and overlaySpws==False and overlayBasebands==False)) |
5239 5372 | print("test2: %s" % (overlayAntennas==False and overlaySpws==False and overlayBasebands==False and (mytime+1==nUniqueTimes or mytime == timerangeList[-1]))) |
5240 5373 | print("test3: %s" % (xant==antennasToPlot[-1] and overlayAntennas==True and overlayTimes==False and overlaySpws==False and overlayBasebands==False)) |
5241 5374 | print("*test4: %s" % ((spwctr==len(spwsToPlot)-1) and (overlaySpws or overlayBasebands) and overlayAntennas==False and overlayTimes==False) ) |
5242 5375 | print(" * = overlaySpws==True" ) |
5243 5376 | print("test5: %s" % (overlayTimes and scanTest)) |
5249 5382 | ((overlaySpws or overlayBasebands) and spwctr+1 >= len(spwsToPlot)) or |
5250 5383 | (overlayAntennas==False and overlaySpws==False and overlayBasebands==False)): |
5251 5384 | mytime += 1 |
5252 5385 | if (debug): |
5253 5386 | print("AT BOTTOM OF LOOP: Incrementing mytime to %d, setting firstUnflaggedAntennaToPlot to 0" % (mytime)) |
5254 5387 | firstUnflaggedAntennaToPlot = 0 # try this |
5255 5388 | doneOverlayTime = False # added on 08-nov-2012 |
5256 5389 | if (overlayBasebands and (uniqueScanNumbers == sorted(scansToPlot))): |
5257 5390 | if (debug): print("Breaking because scans not specified") |
5258 5391 | break |
5259 - | # end of while(mytime) loop |
5392 + | # end of enormous while(mytime) loop endwhile mytime |
5260 5393 | if (redisplay == False): |
5261 5394 | spwctr += 1 |
5262 5395 | if (debug): |
5263 5396 | print( "---------------------------------------- Incrementing spwctr to %d, spwsToPlot=" % (spwctr), spwsToPlot) |
5264 5397 | if (spwctr < len(spwsToPlot)): |
5265 5398 | print("---------------------------------------- ispw = %d" % (spwsToPlot[spwctr])) |
5266 5399 | else: |
5267 5400 | print("---------------------------------------- done the spws in this baseband (%d)" % (baseband)) |
5268 5401 | else: |
5269 5402 | if (debug): |
5383 5516 | else: |
5384 5517 | return(list(msFields).index(token)) |
5385 5518 | |
5386 5519 | def GetFieldNamesForFieldId(u, mymsmd, msFields): |
5387 5520 | if (mymsmd != '' and mymsmd != None): |
5388 5521 | return(mymsmd.namesforfields(u)[0]) |
5389 5522 | else: |
5390 5523 | print("B") |
5391 5524 | return(msFields[u]) |
5392 5525 | |
5393 - | |
5394 5526 | def DrawAntennaNamesForOverlayAntennas(xstartPolLabel, ystartPolLabel, polsToPlot, corr_type, channeldiff, ystartMadLabel, subplotRows, gamp_mad, gamp_std, overlayColors, mysize, ampmarkstyle, markersize, markeredgewidth, msAnt, msFound, antennasToPlot, ampmarkstyle2, xframe, firstFrame, caltableTitle, titlesize, debug=False): |
5395 5527 | if (debug): print("overlayAntennas=True") |
5396 5528 | x0 = xstartPolLabel |
5397 5529 | y0 = ystartPolLabel |
5398 5530 | # draw polarization labels |
5399 5531 | if (debug): print("1) overlayAntennas=True") |
5400 5532 | if (corrTypeToString(corr_type[0]) in polsToPlot): |
5401 5533 | if (channeldiff > 0): |
5402 5534 | pb.text(x0, ystartMadLabel-0.03*subplotRows*0, |
5403 5535 | corrTypeToString(corr_type[0])+' MAD = %.4f, St.Dev = %.4f'%(gamp_mad[0]['mad'],gamp_std[0]['std']), |
5425 5557 | color=overlayColors[0],size=mysize, transform=pb.gca().transAxes) |
5426 5558 | pdesc = pb.plot([x0-0.01], [y0-0.03*subplotRows], '%sk'%ampmarkstyle2, |
5427 5559 | markersize=markersize, scalex=False,scaley=False,markeredgewidth=markeredgewidth) |
5428 5560 | if (debug): print("3) overlayAntennas=True") |
5429 5561 | if (xframe == firstFrame): |
5430 5562 | # draw title including caltable name |
5431 5563 | if (debug): print("4) overlayAntennas=True") |
5432 5564 | pb.text(xstartTitle, ystartTitle, caltableTitle, size=titlesize, color='k', |
5433 5565 | transform=pb.gcf().transFigure) |
5434 5566 | if (debug): print("5) overlayAntennas=True") |
5435 - | DrawAntennaNames(msAnt, antennasToPlot, msFound, mysize) |
5567 + | DrawAntennaNames(msAnt, antennasToPlot, msFound, mysize, overlayColors) |
5436 5568 | if (debug): print("6) overlayAntennas=True") |
5437 5569 | |
5438 5570 | def getTelescopeNameFromCaltable(caltable): |
5439 5571 | mytb = table() |
5440 5572 | mytb.open(caltable) |
5441 5573 | if ('OBSERVATION' in mytb.getkeywords()): |
5442 5574 | observationTable = mytb.getkeyword('OBSERVATION').split()[1] |
5443 5575 | else: |
5444 5576 | observationTable = None |
5445 5577 | mytb.close() |
5702 5834 | bestscan = scans[i] |
5703 5835 | if (verbose): print("bestscan = %s" % (str(bestscan))) |
5704 5836 | mindiff = tdiff |
5705 5837 | if (verbose): |
5706 5838 | print("For timestamp=%.1f, got closest scan = %d, %.0f sec away" %(timestamp, bestscan,mindiff)) |
5707 5839 | if (verbose): print("Calling getWeather()") |
5708 5840 | [conditions,myTimes] = getWeather(vis,bestscan,antenna,verbose,mymsmd) |
5709 5841 | if (verbose): print("Done getWeather()") |
5710 5842 | |
5711 5843 | # convert pressure with unit to the value in mbar |
5844 + | if (verbose): print("conditions = ", conditions) |
5712 5845 | P = myqa.convert(myqa.quantity(conditions['pressure'], conditions['pressure_unit']), 'mbar')['value'] |
5713 - | |
5846 + | if P < 100: ### added 2024Aug |
5847 + | # then the units are wrong, as in early ALMA data had units="Pa" but value was in mbar |
5848 + | P = conditions['pressure'] ### added 2024Aug |
5849 + | casalogPost(verbose,"Ignoring pressure units '%s' because the value implied would be less than 100 mbar."%(conditions['pressure_unit'])) ### added 2024Aug |
5714 5850 | H = conditions['humidity'] |
5715 5851 | T = conditions['temperature']+273.15 |
5716 5852 | if (P <= 0.0): |
5717 5853 | P = 563 |
5718 5854 | if (H <= 0.0): |
5719 5855 | H = 20 |
5720 5856 | else: |
5721 5857 | conditions = {} |
5722 5858 | if (type(mymsmd) != str): |
5723 5859 | if ('elevation' not in list(conditions.keys())): |
5726 5862 | myfieldId = mymsmd.fieldsforscan(bestscan)[0] |
5727 5863 | myscantime = np.mean(mymsmd.timesforscan(bestscan)) |
5728 5864 | mydirection = getRADecForField(vis, myfieldId, verbose) |
5729 5865 | if (verbose): |
5730 5866 | print("myfieldId = %s" % (str(myfieldId))) |
5731 5867 | print("mydirection = %s" % (str(mydirection))) |
5732 5868 | print("Scan = %d, time = %.1f, Field = %d, direction = %s" % (bestscan, myscantime, myfieldId, str(mydirection))) |
5733 5869 | telescopeName = mymsmd.observatorynames()[0] |
5734 5870 | if (len(telescopeName) < 1): |
5735 5871 | telescopeName = 'ALMA' |
5736 - | print("telescope = %s" % (telescopeName)) |
5872 + | if verbose: ### added 2024Aug to prevent lots of spam to the terminal |
5873 + | print("telescope = %s" % (telescopeName)) |
5737 5874 | myazel = computeAzElFromRADecMJD(mydirection, myscantime/86400., telescopeName) |
5738 5875 | conditions['elevation'] = myazel[1] * 180/math.pi |
5739 5876 | conditions['azimuth'] = myazel[0] * 180/math.pi |
5740 5877 | if (verbose): |
5741 5878 | print("Computed elevation = %.1f deg" % (conditions['elevation'])) |
5742 5879 | else: |
5743 - | conditions['elevation'] = 90 |
5880 + | conditions['elevation'] = 45 |
5881 + | casalogPost(verbose,"Using 45 deg elevation since the actual value is unavailable.") |
5744 5882 | bestscan = -1 |
5745 5883 | if (verbose): |
5746 5884 | print("CalcAtm: found elevation=%f (airmass=%.3f) for scan: %s" % (conditions['elevation'],1/np.sin(conditions['elevation']*np.pi/180.), str(bestscan))) |
5747 5885 | print("P,H,T = %f,%f,%f" % (P,H,T)) |
5748 5886 | if (conditions['elevation'] <= 3): |
5749 5887 | print("Using 45 deg elevation instead") |
5750 5888 | airmass = 1.0/math.cos(45*math.pi/180.) |
5751 5889 | else: |
5752 5890 | airmass = 1.0/math.cos((90-conditions['elevation'])*math.pi/180.) |
5753 5891 | |
5754 5892 | # get the elevation of the antenna |
5755 5893 | geodetic_elevation = 5059 |
5756 5894 | if os.path.exists(os.path.join(vis, 'ANTENNA')): |
5757 5895 | with sdutil.table_manager(os.path.join(vis, 'ANTENNA')) as tb: |
5758 5896 | _X, _Y, _Z = (float(i) for i in tb.getcell('POSITION', antenna)) |
5759 5897 | geodetic_elevation = simutil.simutil().xyz2long(_X, _Y, _Z, 'WGS84')[2] |
5760 - | |
5898 + | if verbose: |
5899 + | casalogPost(True,"computed geodetic_elevation from ANTENNA table as %s" % str(geodetic_elevation)) |
5900 + | |
5761 5901 | tropical = 1 |
5762 5902 | midLatitudeSummer = 2 |
5763 5903 | midLatitudeWinter = 3 |
5764 5904 | numchan = len(freqs) |
5765 5905 | # Set the reference freq to be the middle of the middle two channels |
5766 5906 | reffreq=0.5*(freqs[numchan//2-1]+freqs[numchan//2]) |
5767 5907 | originalnumchan = numchan |
5768 5908 | while (numchan > MAX_ATM_CALC_CHANNELS): |
5769 5909 | numchan //= 2 |
5770 5910 | # print("Reducing numchan to ", numchan) |
5771 5911 | chans = list(range(0,originalnumchan,(originalnumchan//numchan))) |
5772 5912 | |
5773 5913 | chansep = (freqs[-1]-freqs[0])/(numchan-1) |
5774 5914 | nbands = 1 |
5775 5915 | if (verbose): print("Opening casac.atmosphere()") |
5776 5916 | myat = atmosphere() |
5777 5917 | if (verbose): print("Opened") |
5778 5918 | fCenter = myqa.quantity(reffreq,'GHz') |
5779 5919 | fResolution = myqa.quantity(chansep,'GHz') |
5780 5920 | fWidth = myqa.quantity(numchan*chansep,'GHz') |
5921 + | h0 = 1.0 # km |
5922 + | dP = 5.0 # mbar |
5923 + | dPm = 1.1 # unitless ratio |
5924 + | # print(H, T, geodetic_elevation,P, midLatitudeWinter, maxAltitude, h0, dP, dPm) |
5781 5925 | myat.initAtmProfile(humidity=H, temperature=myqa.quantity(T, "K"), |
5782 5926 | altitude=myqa.quantity(geodetic_elevation, "m"), |
5783 - | pressure=myqa.quantity(P, 'mbar'), atmType=midLatitudeWinter) |
5927 + | pressure=myqa.quantity(P, 'mbar'), |
5928 + | atmType=midLatitudeWinter, |
5929 + | h0=myqa.quantity(h0,"km"), |
5930 + | maxAltitude=myqa.quantity(maxAltitude,"km"), |
5931 + | dP=myqa.quantity(dP,"mbar"), |
5932 + | dPm=dPm) |
5784 5933 | myat.initSpectralWindow(nbands,fCenter,fWidth,fResolution) |
5785 5934 | myat.setUserWH2O(myqa.quantity(pwvmean,'mm')) |
5786 5935 | |
5787 5936 | # myat.setAirMass() # This does not affect the opacity, but it does effect TebbSky, so do it manually. |
5788 5937 | |
5789 5938 | n = myat.getNumChan() |
5790 5939 | if (verbose): print("numchan = %s" % (str(n))) |
5791 5940 | if ctsys.compare_version('<',[4,0,0]): |
5792 5941 | dry = np.array(myat.getDryOpacitySpec(0)['dryOpacity']) |
5793 5942 | wet = np.array(myat.getWetOpacitySpec(0)['wetOpacity'].value) |
5937 6086 | # otherwise rescale to the automatically-determined range. |
5938 6087 | |
5939 6088 | # chans = 0..N where N=number of channels in the ATM_CALC |
5940 6089 | # channels = 0..X where X=number of channels in the spw, regardless of flagging |
5941 6090 | |
5942 6091 | if (len(chans) != len(channels)): |
5943 6092 | if (chans[1] > chans[0]): |
5944 6093 | atmchanrange = chans[-1]-chans[0] |
5945 6094 | else: |
5946 6095 | atmchanrange = chans[0]-chans[-1] |
5947 - | |
6096 + | if len(channels) == 0: ### added 2024Aug to prevent crash |
6097 + | return(chans) ### added 2024Aug to prevent crash |
5948 6098 | if (channels[1] > channels[0]): |
5949 6099 | chanrange = channels[-1]-channels[0] |
5950 6100 | else: |
5951 6101 | chanrange = channels[0]-channels[-1] |
5952 6102 | |
5953 6103 | newchans = np.array(chans)*chanrange/atmchanrange |
5954 6104 | return(newchans) |
5955 6105 | else: |
5956 6106 | return(chans) |
5957 6107 | |
6050 6200 | |
6051 6201 | def SetNewYLimits(newylimits): |
6052 6202 | # print("Entered SetNewYLimits with ", newylimits ) |
6053 6203 | newrange = newylimits[1]-newylimits[0] |
6054 6204 | if (newrange > 0): |
6055 6205 | pb.ylim([newylimits[0]-0.0*newrange, newylimits[1]+0.0*newrange]) |
6056 6206 | |
6057 6207 | def SetNewXLimits(newxlimits, loc=0): |
6058 6208 | # print("loc=%d: Entered SetNewXLimits with range = %.3f (%f-%f)" % (loc,np.max(newxlimits)-np.min(newxlimits), newxlimits[0], newxlimits[1])) |
6059 6209 | myxrange = np.abs(newxlimits[1]-newxlimits[0]) |
6210 + | if (myxrange == 0): |
6211 + | myxrange = 0.001 |
6060 6212 | mybuffer = 0.01 |
6061 6213 | if (newxlimits[0] < newxlimits[1]): |
6062 6214 | pb.xlim([newxlimits[0]-myxrange*mybuffer,newxlimits[1]+myxrange*mybuffer] ) |
6063 6215 | else: |
6064 6216 | # print("Swapping xlimits order") |
6065 6217 | pb.xlim(newxlimits[1]-myxrange*mybuffer, newxlimits[0]+myxrange*mybuffer) |
6066 6218 | |
6067 6219 | def sloppyMatch(newvalue, mylist, threshold, mytime=None, scansToPlot=[], |
6068 6220 | scansForUniqueTimes=[], myprint=False, whichone=False): |
6069 6221 | """ |
6098 6250 | mymatch = i |
6099 6251 | if (matched == False and myprint==True): |
6100 6252 | print("sloppyMatch: %.0f is not within %.0f of anything in %s" % (newvalue,threshold, str([int(round(b)) for b in mylist]))) |
6101 6253 | elif (myprint==True): |
6102 6254 | print("sloppyMatch: %.0f is within %.0f of something in %s" % (newvalue,threshold, str([int(round(b)) for b in mylist]))) |
6103 6255 | if (whichone == False): |
6104 6256 | return(matched) |
6105 6257 | else: |
6106 6258 | return(matched,mymatch) |
6107 6259 | |
6108 - | def sloppyUnique(t, thresholdSeconds): |
6260 + | # replacement function with extra optional parameter added on 2024Aug |
6261 + | def sloppyUnique(t, thresholdSeconds, cal_scans=None): |
6109 6262 | """ |
6110 - | Takes a list of numbers and returns a list of unique values, subject to a threshold difference. |
6263 + | Takes a list of numbers and returns a list of unique values, subject |
6264 + | to a threshold difference. |
6265 + | cal_scans: if specified, then perform the analysis per scan number |
6111 6266 | """ |
6112 - | # start with the first entry, and only add a new entry if it is more than the threshold from prior |
6113 - | sloppyList = [t[0]] |
6114 - | for i in range(1,len(t)): |
6115 - | keepit = True |
6116 - | for j in range(0,i): |
6117 - | if (abs(t[i]-t[j]) < thresholdSeconds): |
6118 - | keepit = False |
6119 - | if (keepit): |
6120 - | sloppyList.append(t[i]) |
6121 - | # print("sloppyUnique returns %d values from the original %d" % (len(sloppyList), len(t))) |
6267 + | # start with the first entry (t[0]), and only add a new entry if it is more than the threshold from prior |
6268 + | if cal_scans is None: |
6269 + | sloppyList = [t[0]] |
6270 + | for i in range(1,len(t)): |
6271 + | keepit = True |
6272 + | # for j in range(0,i): # prior to PIPE-1519 |
6273 + | # if (abs(t[i]-t[j]) < thresholdSeconds): |
6274 + | for j in range(len(sloppyList)): |
6275 + | if (abs(t[i]-sloppyList[j]) < thresholdSeconds): |
6276 + | keepit = False |
6277 + | if (keepit): |
6278 + | sloppyList.append(t[i]) |
6279 + | else: |
6280 + | # build a list of one time per calibration scan |
6281 + | uniqueCalScans = np.unique(cal_scans) |
6282 + | t = np.array(t) |
6283 + | sloppyList = [] |
6284 + | # print("len(t) = %d, len(cal_scans) = %d" % (len(t),len(cal_scans))) |
6285 + | for k,calscan in enumerate(uniqueCalScans): |
6286 + | times = t[np.where(calscan == cal_scans)] |
6287 + | sloppyList.append(times[0]) |
6288 + | # print "sloppyUnique returns %d values from the original %d" % (len(sloppyList), len(t)) |
6122 6289 | return(sloppyList) |
6123 6290 | |
6291 + | # commented out 2024Aug |
6292 + | #def sloppyUnique(t, thresholdSeconds): |
6293 + | # """ |
6294 + | # Takes a list of numbers and returns a list of unique values, subject to a threshold difference. |
6295 + | # """ |
6296 + | # # start with the first entry, and only add a new entry if it is more than the threshold from prior |
6297 + | # sloppyList = [t[0]] |
6298 + | # for i in range(1,len(t)): |
6299 + | # keepit = True |
6300 + | # for j in range(0,i): |
6301 + | # if (abs(t[i]-t[j]) < thresholdSeconds): |
6302 + | # keepit = False |
6303 + | # if (keepit): |
6304 + | # sloppyList.append(t[i]) |
6305 + | ## print("sloppyUnique returns %d values from the original %d" % (len(sloppyList), len(t))) |
6306 + | # return(sloppyList) |
6307 + | |
6124 6308 | def SetLimits(plotrange, chanrange, newylimits, channels, frequencies, pfrequencies, |
6125 6309 | ampMin, ampMax, xaxis, pxl, chanrangeSetXrange, chanrangePercent=None): |
6126 6310 | """ |
6127 6311 | This is the place where chanrange actually takes effect. |
6128 6312 | """ |
6129 6313 | if (abs(plotrange[0]) > 0 or abs(plotrange[1]) > 0): |
6130 6314 | SetNewXLimits([plotrange[0],plotrange[1]]) |
6131 6315 | if (plotrange[2] == 0 and plotrange[3] == 0): |
6132 6316 | # reset the ylimits based on the channel range shown (selected via plotrange) |
6133 6317 | SetNewYLimits(newylimits) |
6405 6589 | else: |
6406 6590 | # This can remain in axes units since it is the final plot. |
6407 6591 | pb.text(xEdgeLabel, otherEdgeYvalue,'%.0f%%'%(otherEdgeT*100), |
6408 6592 | color=atmcolor, |
6409 6593 | size=mysize, transform=pb.gca().transAxes) |
6410 6594 | pb.text(xEdgeLabel, zeroYValue,'%.0f%%'%(zeroValue*100), |
6411 6595 | color=atmcolor, |
6412 6596 | size=mysize, transform=pb.gca().transAxes) |
6413 6597 | if (lo1 != ''): |
6414 6598 | if (xframe == firstFrame): |
6415 - | pb.text(+0.96-0.08*subplotCols, -0.07*subplotRows, |
6416 - | 'Signal Sideband', color='m', size=mysize, |
6599 + | pb.text(+1.04-0.04*subplotCols, -0.07*subplotRows, |
6600 + | 'Signal SB', color='m', size=mysize, |
6417 6601 | transform=pb.gca().transAxes) |
6418 - | pb.text(-0.08*subplotCols, -0.07*subplotRows, |
6419 - | 'Image Sideband', color='k', size=mysize, |
6602 + | pb.text(-0.03-0.08*subplotCols, -0.07*subplotRows, |
6603 + | 'Image SB', color='k', size=mysize, |
6420 6604 | transform=pb.gca().transAxes) |
6605 + | # pb.text(+0.96-0.08*subplotCols, -0.07*subplotRows, |
6606 + | # 'Signal Sideband', color='m', size=mysize, |
6607 + | # transform=pb.gca().transAxes) |
6608 + | # pb.text(-0.08*subplotCols, -0.07*subplotRows, |
6609 + | # 'Image Sideband', color='k', size=mysize, |
6610 + | # transform=pb.gca().transAxes) |
6421 6611 | return ylim # CAS-8655 |
6422 6612 | |
6423 6613 | def DrawBottomLegendPageCoords(msName, uniqueTimesMytime, mysize, figfile): |
6424 6614 | msName = msName.split('/')[-1] |
6425 6615 | bottomLegend = msName + ' ObsDate=' + utdatestring(uniqueTimesMytime) |
6426 6616 | if (os.path.basename(figfile).find('regression') == 0): |
6427 6617 | regression = True |
6428 6618 | else: |
6429 6619 | regression = False |
6430 6620 | if (regression == False): |