##// END OF EJS Templates
schain
RSS
Specular Meteor module modifications
Julio Valdez -
r840:0301fd9b0cbd
parent child
Show More
Show file before | Show file after | Hide comments
@@ -345,10 +345,11 class ParametersProc(ProcessingUnit):
345 345 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
346 346 noise_timeStep = 4, noise_multiple = 4,
347 347 multDet_timeLimit = 1, multDet_rangeLimit = 3,
348 phaseThresh = 20, SNRThresh = 8,
348 phaseThresh = 20, SNRThresh = 5,
349 349 hmin = 50, hmax=150, azimuth = 0,
350 channel25X = 0, channel20X = 4, channelCentX = 2,
351 channel25Y = 1, channel20Y = 3, channelCentY = 2) :
350 # channel25X = 0, channel20X = 4, channelCentX = 2,
351 # channel25Y = 1, channel20Y = 3, channelCentY = 2,
352 channelPositions = None) :
352 353
353 354 '''
354 355 Function DetectMeteors()
@@ -413,7 +414,14 class ParametersProc(ProcessingUnit):
413 414 Phase0 Phase1 Phase2 Phase3
414 415 TypeError
415 416
416 '''
417 '''
418 #Getting Pairslist
419 if channelPositions == None:
420 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
421 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
422 meteorOps = MeteorOperations()
423 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
424
417 425 #Get Beacon signal
418 426 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
419 427
@@ -421,12 +429,7 class ParametersProc(ProcessingUnit):
421 429 newheis = numpy.where(self.dataOut.heightList>hei_ref)
422 430
423 431 heiRang = self.dataOut.getHeiRange()
424 #Pairs List
425 '''
426 Cambiar este pairslist
427 '''
428 pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
429
432
430 433 # nChannel = self.dataOut.nChannels
431 434 # for i in range(nChannel):
432 435 # if i != centerReceiverIndex:
@@ -467,7 +470,7 class ParametersProc(ProcessingUnit):
467 470 if cohDetection:
468 471 #use coherent detection to get the net power
469 472 cohDet_thresh = cohDet_thresh*numpy.pi/180
470 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
473 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh)
471 474
472 475 #Non-coherent detection!
473 476 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
@@ -498,7 +501,7 class ParametersProc(ProcessingUnit):
498 501 phaseThresh = phaseThresh*numpy.pi/180
499 502 thresh = [phaseThresh, noise_multiple, SNRThresh]
500 503 #Meteor reestimation (Errors N 1, 6, 12, 17)
501 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
504 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
502 505 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
503 506 #Estimation of decay times (Errors N 7, 8, 11)
504 507 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
@@ -507,7 +510,7 class ParametersProc(ProcessingUnit):
507 510 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
508 511 #Calculating Radial Velocity (Error N 15)
509 512 radialStdThresh = 10
510 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
513 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval)
511 514
512 515 if len(listMeteors4) > 0:
513 516 #Setting New Array
@@ -526,10 +529,9 class ParametersProc(ProcessingUnit):
526 529 pairsList.append(pairx)
527 530 pairsList.append(pairy)
528 531
529 meteorOps = MeteorOperations()
530 532 jph = numpy.array([0,0,0,0])
531 533 h = (hmin,hmax)
532 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
534 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
533 535
534 536 # #Calculate AOA (Error N 3, 4)
535 537 # #JONES ET AL. 1998
@@ -552,10 +554,12 class ParametersProc(ProcessingUnit):
552 554
553 555 if arrayParameters == None:
554 556 self.dataOut.flagNoData = True
557 else:
558 self.dataOut.flagNoData = True
555 559
556 560 return
557 561
558 def CorrectMeteorPhases(self, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45):
562 def CorrectMeteorPhases(self, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
559 563
560 564 arrayParameters = self.dataOut.data_param
561 565 pairsList = []
@@ -566,12 +570,19 class ParametersProc(ProcessingUnit):
566 570 jph = numpy.zeros(4)
567 571
568 572 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
569 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
573 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
574 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
570 575
571 576 meteorOps = MeteorOperations()
577 if channelPositions == None:
578 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
579 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
580
581 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
572 582 h = (hmin,hmax)
583
584 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
573 585
574 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
575 586 self.dataOut.data_param = arrayParameters
576 587 return
577 588
@@ -641,10 +652,11 class ParametersProc(ProcessingUnit):
641 652 phaseArrival = phaseInt.reshape(phaseInt.size)
642 653
643 654 #Dealias
644 indAlias = numpy.where(phaseArrival > numpy.pi)
645 phaseArrival[indAlias] -= 2*numpy.pi
646 indAlias = numpy.where(phaseArrival < -numpy.pi)
647 phaseArrival[indAlias] += 2*numpy.pi
655 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
656 # indAlias = numpy.where(phaseArrival > numpy.pi)
657 # phaseArrival[indAlias] -= 2*numpy.pi
658 # indAlias = numpy.where(phaseArrival < -numpy.pi)
659 # phaseArrival[indAlias] += 2*numpy.pi
648 660
649 661 return phaseDiff, phaseArrival
650 662
@@ -2293,17 +2305,25 class PhaseCalibration(Operation):
2293 2305
2294 2306 return False
2295 2307
2296 def __getGammas(self, pairs, k, d, phases):
2308 def __getGammas(self, pairs, d, phases):
2297 2309 gammas = numpy.zeros(2)
2298 2310
2299 2311 for i in range(len(pairs)):
2300 2312
2301 2313 pairi = pairs[i]
2302 2314
2315 phip3 = phases[:,pairi[1]]
2316 d3 = d[pairi[1]]
2317 phip2 = phases[:,pairi[0]]
2318 d2 = d[pairi[0]]
2303 2319 #Calculating gamma
2304 jdcos = phases[:,pairi[1]]/(k*d[pairi[1]])
2305 jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]])))
2306
2320 # jdcos = alp1/(k*d1)
2321 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2322 jgamma = -phip2*d3/d2 - phip3
2323 jgamma = numpy.angle(numpy.exp(1j*jgamma))
2324 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
2325 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
2326
2307 2327 #Revised distribution
2308 2328 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
2309 2329
@@ -2346,7 +2366,7 class PhaseCalibration(Operation):
2346 2366 def __gauss_function(self, t, coeffs):
2347 2367
2348 2368 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2349
2369
2350 2370 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2351 2371 meteorOps = MeteorOperations()
2352 2372 nchan = 4
@@ -2379,14 +2399,14 class PhaseCalibration(Operation):
2379 2399 for iy in range(int(nstepsy)):
2380 2400 for ix in range(int(nstepsx)):
2381 2401 jph[pairy[1]] = alpha_y[iy]
2382 jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]]
2402 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2383 2403
2384 2404 jph[pairx[1]] = alpha_x[ix]
2385 jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]]
2405 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2386 2406
2387 2407 jph_array[:,ix,iy] = jph
2388 2408
2389 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph)
2409 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2390 2410 error = meteorsArray1[:,-1]
2391 2411 ind1 = numpy.where(error==0)[0]
2392 2412 penalty[ix,iy] = ind1.size
@@ -2402,7 +2422,7 class PhaseCalibration(Operation):
2402 2422 return phOffset
2403 2423
2404 2424
2405 def run(self, dataOut, hmin, hmax, direction25X=-1, direction20X=1, direction25Y=-1, direction20Y=1, nHours = 1):
2425 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
2406 2426
2407 2427 dataOut.flagNoData = True
2408 2428 self.__dataReady = False
@@ -2435,7 +2455,14 class PhaseCalibration(Operation):
2435 2455 azimuth = 0
2436 2456 h = (hmin, hmax)
2437 2457 pairs = ((0,1),(2,3))
2438 distances = [direction25X*2.5*lamb, direction20X*2*lamb, direction25Y*2.5*lamb, direction20Y*2*lamb]
2458
2459 if channelPositions == None:
2460 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2461 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2462 meteorOps = MeteorOperations()
2463 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2464
2465 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2439 2466
2440 2467 meteorsArray = self.__buffer
2441 2468 error = meteorsArray[:,-1]
@@ -2446,7 +2473,7 class PhaseCalibration(Operation):
2446 2473 phases = meteorsArray[:,8:12]
2447 2474
2448 2475 #Calculate Gammas
2449 gammas = self.__getGammas(pairs, k, distances, phases)
2476 gammas = self.__getGammas(pairs, distances, phases)
2450 2477 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2451 2478 #Calculate Phases
2452 2479 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
@@ -2464,7 +2491,7 class MeteorOperations():
2464 2491
2465 2492 return
2466 2493
2467 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph):
2494 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
2468 2495
2469 2496 arrayParameters = arrayParameters0.copy()
2470 2497 hmin = h[0]
@@ -2476,7 +2503,7 class MeteorOperations():
2476 2503 error = arrayParameters[:,-1]
2477 2504 phases = -arrayParameters[:,8:12] + jph
2478 2505 # phases = numpy.unwrap(phases)
2479 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2506 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
2480 2507
2481 2508 #Calculate Heights (Error N 13 and 14)
2482 2509 error = arrayParameters[:,-1]
@@ -2491,10 +2518,10 class MeteorOperations():
2491 2518
2492 2519 return arrayParameters
2493 2520
2494 def getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2521 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
2495 2522
2496 2523 arrayAOA = numpy.zeros((phases.shape[0],3))
2497 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2524 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
2498 2525
2499 2526 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2500 2527 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
@@ -2514,7 +2541,7 class MeteorOperations():
2514 2541 error[indInvalid] = 4
2515 2542 return arrayAOA, error
2516 2543
2517 def __getDirectionCosines(self, arrayPhase, pairsList):
2544 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2518 2545
2519 2546 #Initializing some variables
2520 2547 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
@@ -2525,21 +2552,23 class MeteorOperations():
2525 2552
2526 2553
2527 2554 for i in range(2):
2555 ph0 = arrayPhase[:,pairsList[i][0]]
2556 ph1 = arrayPhase[:,pairsList[i][1]]
2557 d0 = distances[pairsList[i][0]]
2558 d1 = distances[pairsList[i][1]]
2559
2560 ph0_aux = ph0 + ph1
2561 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
2562 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
2563 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
2528 2564 #First Estimation
2529 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2530 #Dealias
2531 indcsi = numpy.where(phi0_aux > numpy.pi)
2532 phi0_aux[indcsi] -= 2*numpy.pi
2533 indcsi = numpy.where(phi0_aux < -numpy.pi)
2534 phi0_aux[indcsi] += 2*numpy.pi
2535 #Direction Cosine 0
2536 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2565 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
2537 2566
2538 2567 #Most-Accurate Second Estimation
2539 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2568 phi1_aux = ph0 - ph1
2540 2569 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2541 2570 #Direction Cosine 1
2542 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2571 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
2543 2572
2544 2573 #Searching the correct Direction Cosine
2545 2574 cosdir0_aux = cosdir0[:,i]
@@ -2557,7 +2586,7 class MeteorOperations():
2557 2586 cosdirY = cosdir[:,1]
2558 2587
2559 2588 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2560 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2589 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
2561 2590 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2562 2591
2563 2592 return angles
@@ -2596,4 +2625,88 class MeteorOperations():
2596 2625 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2597 2626 error[indInvalid1] = 13
2598 2627
2599 return heights, error No newline at end of file
2628 return heights, error
2629
2630 def getPhasePairs(self, channelPositions):
2631 chanPos = numpy.array(channelPositions)
2632 listOper = list(itertools.combinations(range(5),2))
2633
2634 distances = numpy.zeros(4)
2635 axisX = []
2636 axisY = []
2637 distX = numpy.zeros(3)
2638 distY = numpy.zeros(3)
2639 ix = 0
2640 iy = 0
2641
2642 pairX = numpy.zeros((2,2))
2643 pairY = numpy.zeros((2,2))
2644
2645 for i in range(len(listOper)):
2646 pairi = listOper[i]
2647
2648 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
2649
2650 if posDif[0] == 0:
2651 axisY.append(pairi)
2652 distY[iy] = posDif[1]
2653 iy += 1
2654 elif posDif[1] == 0:
2655 axisX.append(pairi)
2656 distX[ix] = posDif[0]
2657 ix += 1
2658
2659 for i in range(2):
2660 if i==0:
2661 dist0 = distX
2662 axis0 = axisX
2663 else:
2664 dist0 = distY
2665 axis0 = axisY
2666
2667 side = numpy.argsort(dist0)[:-1]
2668 axis0 = numpy.array(axis0)[side,:]
2669 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
2670 axis1 = numpy.unique(numpy.reshape(axis0,4))
2671 side = axis1[axis1 != chanC]
2672 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
2673 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
2674 if diff1<0:
2675 chan2 = side[0]
2676 d2 = numpy.abs(diff1)
2677 chan1 = side[1]
2678 d1 = numpy.abs(diff2)
2679 else:
2680 chan2 = side[1]
2681 d2 = numpy.abs(diff2)
2682 chan1 = side[0]
2683 d1 = numpy.abs(diff1)
2684
2685 if i==0:
2686 chanCX = chanC
2687 chan1X = chan1
2688 chan2X = chan2
2689 distances[0:2] = numpy.array([d1,d2])
2690 else:
2691 chanCY = chanC
2692 chan1Y = chan1
2693 chan2Y = chan2
2694 distances[2:4] = numpy.array([d1,d2])
2695 # axisXsides = numpy.reshape(axisX[ix,:],4)
2696 #
2697 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
2698 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
2699 #
2700 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
2701 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
2702 # channel25X = int(pairX[0,ind25X])
2703 # channel20X = int(pairX[1,ind20X])
2704 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
2705 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
2706 # channel25Y = int(pairY[0,ind25Y])
2707 # channel20Y = int(pairY[1,ind20Y])
2708
2709 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
2710 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
2711
2712 return pairslist, distances No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now