casa5
Public
Commits
Ville Suoranta authored 10425691988 Merge
| 1097 1097 | |
| 1098 1098 | |
| 1099 1099 | |
| 1100 1100 | def sensitivity(self, freq, bandwidth, etime, elevation, pwv=None, |
| 1101 1101 | telescope=None, diam=None, nant=None, |
| 1102 1102 | antennalist=None, |
| 1103 1103 | doimnoise=None, |
| 1104 1104 | integration=None,debug=None, |
| 1105 1105 | method="tsys-atm",tau0=None,t_sky=None): |
| 1106 1106 | |
| 1107 + | if qa.convert(elevation,"deg")['value']<3: |
| 1108 + | self.msg("Elevation < ALMA limit of 3 deg",priority="error") |
| 1109 + | return False |
| 1110 + | |
| 1107 1111 | import glob |
| 1108 1112 | tmpname="tmp"+str(os.getpid()) |
| 1109 1113 | i=0 |
| 1110 1114 | while i<10 and len(glob.glob(tmpname+"*"))>0: |
| 1111 1115 | tmpname="tmp"+str(os.getpid())+str(i) |
| 1112 1116 | i=i+1 |
| 1113 1117 | if i>=9: |
| 1114 1118 | xx=glob.glob(tmpname+"*") |
| 1115 1119 | for k in range(len(xx)): |
| 1116 1120 | if os.path.isdir(xx[k]): |
| 1159 1163 | resl=qa.convert(tail,"arcsec")['value'] |
| 1160 1164 | repodir=os.getenv("CASAPATH").split(' ')[0]+"/data/alma/simmos/" |
| 1161 1165 | if os.path.exists(repodir): |
| 1162 1166 | confnum=(2.867-pl.log10(resl*1000*qa.convert(freq,"GHz")['value']/672.))/0.0721 |
| 1163 1167 | confnum=max(1,min(28,confnum)) |
| 1164 1168 | conf=str(int(round(confnum))) |
| 1165 1169 | if len(conf)<2: conf='0'+conf |
| 1166 1170 | antennalist=repodir+"alma.out"+conf+".cfg" |
| 1167 1171 | self.msg("converted resolution to antennalist "+antennalist) |
| 1168 1172 | |
| 1173 + | repodir = os.getenv("CASAPATH").split(' ')[0] + "/data/alma/simmos/" |
| 1174 + | if not os.path.exists(antennalist) and \ |
| 1175 + | os.path.exists(repodir+antennalist): |
| 1176 + | antennalist = repodir + antennalist |
| 1169 1177 | if os.path.exists(antennalist): |
| 1170 1178 | stnx, stny, stnz, stnd, padnames, telescope, posobs = self.readantenna(antennalist) |
| 1171 1179 | else: |
| 1172 1180 | self.msg("antennalist "+antennalist+" not found",priority="error") |
| 1173 1181 | return False |
| 1174 1182 | |
| 1183 + | nant=len(stnx) |
| 1175 1184 | # diam is only used as a test below, not quantitatively |
| 1176 1185 | diam = pl.average(stnd) |
| 1177 1186 | antnames=padnames |
| 1178 1187 | |
| 1179 1188 | |
| 1180 1189 | if (telescope==None or diam==None): |
| 1181 1190 | self.msg("Telescope name or antenna diameter have not been set.",priority="error") |
| 1182 1191 | return False |
| 1183 1192 | |
| 1184 1193 | # copied from task_simdata: |
| 1192 1201 | # start is center of first channel. for nch=1, that equals center |
| 1193 1202 | model_start=qa.quantity(freq) |
| 1194 1203 | |
| 1195 1204 | stokes, feeds = self.polsettings(telescope) |
| 1196 1205 | sm.setspwindow(spwname="band1", freq=qa.tos(model_start), |
| 1197 1206 | deltafreq=qa.tos(model_width), |
| 1198 1207 | freqresolution=qa.tos(model_width), |
| 1199 1208 | nchannels=model_nchan, stokes=stokes) |
| 1200 1209 | sm.setfeed(mode=feeds, pol=['']) |
| 1201 1210 | |
| 1202 - | sm.setlimits(shadowlimit=0.01, elevationlimit='10deg') |
| 1211 + | sm.setlimits(shadowlimit=0.01, elevationlimit='3deg') |
| 1203 1212 | sm.setauto(0.0) |
| 1204 1213 | |
| 1205 - | obslat=qa.convert(obs['m1'],'deg') |
| 1214 + | obslat=qa.convert(posobs['m1'],'deg') |
| 1206 1215 | dec=qa.add(obslat, qa.add(qa.quantity("90deg"),qa.mul(elevation,-1))) |
| 1207 1216 | |
| 1208 1217 | sm.setfield(sourcename="src1", |
| 1209 1218 | sourcedirection="J2000 00:00:00.00 "+qa.angle(dec)[0], |
| 1210 1219 | calcode="OBJ", distance='0m') |
| 1211 1220 | reftime = me.epoch('TAI', "2012/01/01/00:00:00") |
| 1212 1221 | if integration==None: |
| 1213 1222 | integration=qa.mul(etime,0.01) |
| 1214 1223 | self.msg("observing for "+qa.tos(etime)+" with integration="+qa.tos(integration)) |
| 1215 1224 | sm.settimes(integrationtime=integration, usehourangle=True, |
| 1281 1290 | imnoise=(stats["rms"][0]) |
| 1282 1291 | else: |
| 1283 1292 | imnoise=0. |
| 1284 1293 | |
| 1285 1294 | nint = qa.convert(etime,'s')['value'] / qa.convert(integration,'s')['value'] |
| 1286 1295 | nbase= 0.5*nant*(nant-1) |
| 1287 1296 | |
| 1288 1297 | if os.path.exists(tmpname+".T.cal"): |
| 1289 1298 | tb.open(tmpname+".T.cal") |
| 1290 1299 | gain=tb.getcol("CPARAM") |
| 1300 + | flag=tb.getcol("FLAG") |
| 1301 + | # RI TODO average instead of first? |
| 1291 1302 | tb.done() |
| 1292 1303 | # gain is per ANT so square for per baseline; |
| 1293 1304 | # pick a gain from about the middle of the track |
| 1294 - | # TODO average over the track instead? |
| 1295 - | noiseperbase=1./(gain[0][0][0.5*nint*nant].real)**2 |
| 1305 + | # needs to be unflagged! |
| 1306 + | z=pl.where(flag[0][0]==False)[0] |
| 1307 + | nunflagged=len(z) |
| 1308 + | # noiseperbase=1./(gain[0][0][0.5*nint*nant].real)**2 |
| 1309 + | noiseperbase=1./(gain[0][0][z[nunflagged/2]].real)**2 |
| 1296 1310 | else: |
| 1297 1311 | noiseperbase=0. |
| 1298 1312 | |
| 1299 1313 | theoreticalnoise=noiseperbase/pl.sqrt(nint*nbase*2) # assume 2-poln |
| 1300 - | |
| 1301 - | if debug==None: |
| 1314 + | |
| 1315 + | if debug!=True: |
| 1302 1316 | xx=glob.glob(tmpname+"*") |
| 1303 1317 | for k in range(len(xx)): |
| 1304 1318 | if os.path.isdir(xx[k]): |
| 1305 1319 | cu.removetable(xx[k]) |
| 1306 1320 | else: |
| 1307 1321 | os.remove(xx[k]) |
| 1308 1322 | |
| 1309 1323 | if doimnoise: |
| 1310 1324 | return theoreticalnoise , imnoise |
| 1311 1325 | else: |