Commits
2060 2060 | Int os; |
2061 2061 | if (!aTerm.isNoOp()) os=aTerm.getOversampling(); |
2062 2062 | else if (!wTerm.isNoOp()) os=wTerm.getOversampling(); |
2063 2063 | else os=psTerm.getOversampling(); |
2064 2064 | return os; |
2065 2065 | } |
2066 2066 | void AWConvFunc::makeAConvFunc(Array<Complex>& convFunc, |
2067 2067 | Array<Complex>& wtConvFunc, CoordinateSystem& csys, |
2068 2068 | Vector<Int>& asupport, Int& npix, |
2069 2069 | const Vector<Double>& freqlist, const Bool doSquint, |
2070 - | const Double& pa){ |
2070 + | const Double& pa, const bool isSingleField){ |
2071 2071 | |
2072 2072 | //Assuming first freq is lowest |
2073 2073 | Quantity cell; |
2074 2074 | Int convnx=0; |
2075 2075 | if(freqlist.nelements() ==0) |
2076 2076 | throw(AipsError("Programmer error: No frequencies has been sent for A-Terms construnction")); |
2077 2077 | asupport.resize(freqlist.nelements()); |
2078 2078 | if(!atermMaker_p) |
2079 2079 | throw(AipsError("Programmer Error: AWConvFunc has to be constructed with a EVLAAperture")); |
2080 2080 | String bandname = EVLAAperture::getVLABandName( |
2081 2081 | freqlist[int(freqlist.nelements() / 2)], |
2082 2082 | atermMaker_p->getTelescopeName()); |
2083 2083 | // cerr << "BANDNAME " << bandname << " telescip " << |
2084 2084 | //atermMaker_p->getTelescopeName() << " csys tel " << |
2085 2085 | //csys.obsInfo().telescope() << endl; |
2086 2086 | std::tie(cell,convnx)=getBeamCellSize(bandname); |
2087 - | //widen it to twice first sidelobe |
2088 - | convnx *=2; |
2089 2087 | convnx = Int(ceil(Float(convnx) / 8.0)) * 8; |
2090 2088 | //////////////////// |
2091 2089 | //cerr << "@@@cell " << cell << " pbnpix " << convnx << " imnpix "<< npix << " imcell"<< csys_p.increment()<< endl; |
2092 2090 | csys_p=csys; |
2093 2091 | Vector<String> units=csys_p.worldAxisUnits(); |
2094 2092 | Vector<Double> incr=csys_p.increment(); |
2095 2093 | Double inpFov=fabs(incr[0]*npix); |
2096 - | |
2097 - | |
2098 - | incr[0]=cell.get(units[0]).getValue(); |
2099 - | incr[1]=cell.get(units[1]).getValue(); |
2100 - | |
2101 - | //cerr << "###inp fov" << inpFov << " pbfov " << fabs(incr[0]*Double(convnx)) << endl; |
2102 - | Double pbFov= fabs(incr[0]*Double(convnx)); |
2103 - | /*if(inpFov > (pbFov/8.0) && inpFov< pbFov/2.0 ){ |
2104 - | convnx = Int(floor(Double(convnx) * inpFov / pbFov / 8.0)) * 8; |
2105 - | pbFov = fabs(incr[0] * Double(convnx)); |
2106 - | cerr << "pbFov " << pbFov << " inpFov " << inpFov << " convnx " |
2107 - | << convnx << endl; |
2094 + | Double pbFov= fabs(cell.get(units[0]).getValue()*Double(convnx)); |
2095 + | cerr << "@@inpfov " << inpFov << " pbFov " << pbFov << endl; |
2096 + | if (inpFov > 0.125 * pbFov) { |
2097 + | incr[0] = cell.get(units[0]).getValue(); |
2098 + | incr[1] = cell.get(units[1]).getValue(); |
2099 + | //For small fov post fft resampling is not good enough...create alarge fov |
2100 + | // here to resample finer at calculation |
2101 + | if(inpFov < pbFov || isSingleField){ |
2102 + | incr[0]*=2.0; |
2103 + | incr[1]*=2.0; |
2104 + | convnx *=2.0; |
2105 + | pbFov *=4.0; |
2106 + | } |
2107 + | npix = convnx; // return that npix |
2108 + | |
2109 + | } else { |
2110 + | // Very small image inside mainlobes |
2111 + | // use the image incr rather than the minimum required |
2112 + | // have to keep convnx bigger as it BeamCalc is unhappy to generate |
2113 + | // beam for small fields no need to calculate into the lobes thus the |
2114 + | // factor 2 |
2115 + | convnx = Int(ceil(fabs(Float(convnx) / incr[0] * |
2116 + | cell.get(units[0]).getValue()) / |
2117 + | 16.0)) * |
2118 + | 8; |
2119 + | npix = int(std::ceil(inpFov / pbFov * Double(convnx) / 2.0)) * 2; |
2120 + | pbFov = fabs(incr[0]) * convnx; |
2121 + | // cerr << "$$$$ npix " << npix << " cnx " << convnx << endl; |
2108 2122 | } |
2109 - | */ |
2110 - | |
2111 - | npix=convnx; |
2112 - | //if((inpFov/pbFov) < 0.5){ |
2113 - | // npix=int(std::ceil(inpFov/pbFov*Double(convnx)/2.0))*2; |
2114 - | // cerr << "$$$$ npix " << npix << " cnx " << convnx << endl; |
2115 - | //} |
2116 - | Vector<Int> stoks={Stokes::RR, Stokes::RL, Stokes::LR, Stokes::LL}; |
2123 + | |
2124 + | Vector<Int> stoks={Stokes::RR, Stokes::RL, Stokes::LR, Stokes::LL}; |
2117 2125 | StokesCoordinate stokesCoords(stoks); |
2118 2126 | Quantum<Vector<Double> > freqs(freqlist, "Hz"); |
2119 2127 | SpectralCoordinate specCoord(MFrequency::TOPO, freqs); |
2120 2128 | CoordinateSystem csysA=csys_p; |
2121 2129 | csysA.setIncrement(incr); |
2122 2130 | Vector<Double> refpix=csysA.referencePixel(); |
2123 2131 | refpix[0]=refpix[1]=Double(convnx)/2.0; |
2124 2132 | csysA.setReferencePixel(refpix); |
2125 2133 | csysA.replaceCoordinate(stokesCoords, 1); |
2126 2134 | csysA.replaceCoordinate(specCoord,2); |
2181 2189 | cerr << "@@@neArr shape "<< newArr.shape() << endl; |
2182 2190 | m.putMiddle(arr, newArr); |
2183 2191 | newArr.resize(); |
2184 2192 | newArr=m.resampleViaFFT(wtArr, fac, fac); |
2185 2193 | wtArr.set(0.0); |
2186 2194 | m.putMiddle(wtArr, newArr); |
2187 2195 | |
2188 2196 | |
2189 2197 | }*/ |
2190 2198 | //cerr << "Post FT MAX arr "<< max(wtArr) << " min "<< min(wtArr) << endl; |
2191 - | |
2192 - | supportAndNormalizeAFunc(support, arr, wtArr); |
2199 + | cerr << "inpFov/pbFov " << (inpFov / pbFov) << endl; |
2200 + | supportAndNormalizeAFunc(support, arr, wtArr, (inpFov/pbFov >1) ); |
2193 2201 | //cerr << "Post Norm MAX arr "<< max(arr) << " min "<< min(arr) << endl; |
2194 2202 | if(k==0){ |
2195 2203 | |
2196 2204 | IPosition tmpShape=shp; |
2197 2205 | tmpShape[3]=freqlist.nelements(); |
2198 2206 | //This is going to be returned to allow resacling etc later |
2199 2207 | //csys=refim::SynthesisUtils::makeUVCoords(csysA, tmpShape); |
2200 2208 | IPosition convShp(4, 2*support, 2*support, 4, freqlist.nelements()); |
2201 2209 | refpix[0]=refpix[1]=support; |
2202 2210 | csys.setReferencePixel(refpix); |
2216 2224 | |
2217 2225 | } |
2218 2226 | //cerr << "Post Slicing MAX arr "<< max(wtConvFunc) << " min "<< min(wtConvFunc) << endl; |
2219 2227 | //cerr << "@@@support " << max(asupport) << endl; |
2220 2228 | } |
2221 2229 | |
2222 2230 | void AWConvFunc::makeAWConvFunc(Array<Complex>& convFunc, |
2223 2231 | Array<Complex>& wtconv, |
2224 2232 | CoordinateSystem& csys, |
2225 2233 | Matrix<Int>& awsupport, Int& npix, const Vector<Double>& freqlist, |
2226 - | const Vector<Double>& wVals, const Bool dosquint, const Double& pa){ |
2234 + | const Vector<Double>& wVals, const Bool dosquint, const Double& pa, |
2235 + | const bool isSingleField){ |
2227 2236 | |
2228 2237 | |
2229 2238 | Array<Complex> pbFT; |
2230 2239 | Array<Complex> pbWFT; |
2231 2240 | Vector<Int> aTsup; |
2232 - | makeAConvFunc(pbFT, pbWFT, csys, aTsup, npix, freqlist, dosquint, pa); |
2241 + | makeAConvFunc(pbFT, pbWFT, csys, aTsup, npix, freqlist, dosquint, pa, isSingleField); |
2233 2242 | //cerr << "In AW wtConv "<< max(pbWFT) << " min "<< min(pbWFT) << endl; |
2234 2243 | /////TESTOO |
2235 2244 | //{ |
2236 2245 | // cerr << "ATSUP " << aTsup << endl; |
2237 2246 | //PagedImage<Complex> booltaq(pbFT.shape(), csys, "BOOLTAQ"); |
2238 2247 | //booltaq.put(pbFT); |
2239 2248 | //} |
2240 2249 | ///// |
2241 2250 | WPConvFunc wpc; |
2242 2251 | Cube<Complex> wCon; |
2408 2417 | it2.array() *= mult; |
2409 2418 | //cerr << std::setprecision(12) << "it2.shape " << it2.array().shape() << " //pos " << it2.pos() << "chan " << chan << " sumUnder " << sumUnder(chan) << " w " << w << endl; |
2410 2419 | } |
2411 2420 | it2.next(); |
2412 2421 | } |
2413 2422 | |
2414 2423 | } |
2415 2424 | return True; |
2416 2425 | |
2417 2426 | } |
2418 - | Bool AWConvFunc::supportAndNormalizeAFunc(Int& sup, Array<Complex>& conv, Array<Complex>& wtconv){ |
2427 + | Bool AWConvFunc::supportAndNormalizeAFunc(Int& sup, Array<Complex>& conv, Array<Complex>& wtconv, const bool isLarge){ |
2419 2428 | sup=-1; |
2420 2429 | IPosition begin(4, 0, 0, 0, 0); |
2421 2430 | IPosition end=conv.shape()-1; |
2422 2431 | Int convSize=conv.shape()[0]; |
2423 2432 | end[2]=0; |
2424 2433 | Matrix<Complex> convPlane(wtconv(begin,end).reform(IPosition(2, convSize, convSize))); |
2425 2434 | Float maxAbsConvFunc, minAbsConvFunc; |
2426 2435 | IPosition minpos, maxpos; |
2427 2436 | minMax(minAbsConvFunc, maxAbsConvFunc, minpos, maxpos, amplitude(convPlane)); |
2428 2437 | Bool found=false; |
2429 2438 | Int trial=0; |
2430 - | for (trial=convSize/2-2; trial>0; trial--) { |
2431 - | //Searching down a diagonal |
2432 - | if(abs(convPlane(convSize/2-trial,convSize/2-trial)) > (5e-3*maxAbsConvFunc) ) { |
2433 - | found=true; |
2434 - | trial=Int(sqrt(2.0*Float(trial*trial))); |
2435 - | |
2436 - | break; |
2439 + | Float suplimit = 5e-3; |
2440 + | if (isLarge) |
2441 + | suplimit = 1e-2; |
2442 + | cerr << "###SUPLIMIT " << suplimit << endl; |
2443 + | for (trial = convSize / 2 - 2; trial > 0; trial--) { |
2444 + | // Searching down a diagonal |
2445 + | if (abs(convPlane(convSize / 2 - trial, convSize / 2 - trial)) > |
2446 + | (suplimit * maxAbsConvFunc)) { |
2447 + | found = true; |
2448 + | trial = Int(sqrt(2.0 * Float(trial * trial))); |
2449 + | |
2450 + | break; |
2451 + | } |
2437 2452 | } |
2438 - | } |
2439 2453 | if(!found) { |
2440 - | if((maxAbsConvFunc-minAbsConvFunc) > (5e-3*maxAbsConvFunc)) |
2454 + | if((maxAbsConvFunc-minAbsConvFunc) > (suplimit*maxAbsConvFunc)) |
2441 2455 | found=true; |
2442 2456 | // if it drops by more than 2 magnitudes per pixel |
2443 2457 | trial= (convSize >10) ? 5 : (convSize/2 - 4); |
2444 2458 | } |
2445 2459 | |
2446 2460 | |
2447 2461 | if(found) { |
2448 2462 | if(trial < 5) |
2449 2463 | trial= ( (convSize >10 ) ? 5 : (convSize/2 - 4)); |
2450 2464 | sup=trial+1; |
2477 2491 | pbIt.next(); |
2478 2492 | wtIt.next(); |
2479 2493 | |
2480 2494 | } |
2481 2495 | } |
2482 2496 | |
2483 2497 | return found; |
2484 2498 | } |
2485 2499 | std::pair<Quantity, int> AWConvFunc::getBeamCellSize(const String& band){ |
2486 2500 | //testoo |
2487 - | Quantity fov(0.048,"rad"); //fov at 1 GHz for VLA |
2501 + | Quantity fov(0.024,"rad"); //fov at 1 GHz for VLA |
2488 2502 | Quantity cell=fov/256; |
2489 2503 | if(band=="EVLA_S") |
2490 2504 | cell=cell/2.0; |
2491 2505 | else if(band=="EVLA_C" || band=="VLA_C") |
2492 2506 | cell=cell/4.0; |
2493 2507 | else if(band=="EVLA_X" || band=="VLA_X") |
2494 2508 | cell=cell/8.0; |
2495 2509 | else if(band=="EVLA_U" || band=="VLA_U") |
2496 2510 | cell=cell/12.0; |
2497 2511 | else if(band=="EVLA_K" || band == "VLA_K") |