##// 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 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
345 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
346 noise_timeStep = 4, noise_multiple = 4,
346 noise_timeStep = 4, noise_multiple = 4,
347 multDet_timeLimit = 1, multDet_rangeLimit = 3,
347 multDet_timeLimit = 1, multDet_rangeLimit = 3,
348 phaseThresh = 20, SNRThresh = 8,
348 phaseThresh = 20, SNRThresh = 5,
349 hmin = 50, hmax=150, azimuth = 0,
349 hmin = 50, hmax=150, azimuth = 0,
350 channel25X = 0, channel20X = 4, channelCentX = 2,
350 # channel25X = 0, channel20X = 4, channelCentX = 2,
351 channel25Y = 1, channel20Y = 3, channelCentY = 2) :
351 # channel25Y = 1, channel20Y = 3, channelCentY = 2,
352 channelPositions = None) :
352
353
353 '''
354 '''
354 Function DetectMeteors()
355 Function DetectMeteors()
@@ -413,7 +414,14 class ParametersProc(ProcessingUnit):
413 Phase0 Phase1 Phase2 Phase3
414 Phase0 Phase1 Phase2 Phase3
414 TypeError
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 #Get Beacon signal
425 #Get Beacon signal
418 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
426 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
419
427
@@ -421,12 +429,7 class ParametersProc(ProcessingUnit):
421 newheis = numpy.where(self.dataOut.heightList>hei_ref)
429 newheis = numpy.where(self.dataOut.heightList>hei_ref)
422
430
423 heiRang = self.dataOut.getHeiRange()
431 heiRang = self.dataOut.getHeiRange()
424 #Pairs List
432
425 '''
426 Cambiar este pairslist
427 '''
428 pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
429
430 # nChannel = self.dataOut.nChannels
433 # nChannel = self.dataOut.nChannels
431 # for i in range(nChannel):
434 # for i in range(nChannel):
432 # if i != centerReceiverIndex:
435 # if i != centerReceiverIndex:
@@ -467,7 +470,7 class ParametersProc(ProcessingUnit):
467 if cohDetection:
470 if cohDetection:
468 #use coherent detection to get the net power
471 #use coherent detection to get the net power
469 cohDet_thresh = cohDet_thresh*numpy.pi/180
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 #Non-coherent detection!
475 #Non-coherent detection!
473 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
476 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
@@ -498,7 +501,7 class ParametersProc(ProcessingUnit):
498 phaseThresh = phaseThresh*numpy.pi/180
501 phaseThresh = phaseThresh*numpy.pi/180
499 thresh = [phaseThresh, noise_multiple, SNRThresh]
502 thresh = [phaseThresh, noise_multiple, SNRThresh]
500 #Meteor reestimation (Errors N 1, 6, 12, 17)
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 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
505 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
503 #Estimation of decay times (Errors N 7, 8, 11)
506 #Estimation of decay times (Errors N 7, 8, 11)
504 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
507 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
@@ -507,7 +510,7 class ParametersProc(ProcessingUnit):
507 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
510 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
508 #Calculating Radial Velocity (Error N 15)
511 #Calculating Radial Velocity (Error N 15)
509 radialStdThresh = 10
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 if len(listMeteors4) > 0:
515 if len(listMeteors4) > 0:
513 #Setting New Array
516 #Setting New Array
@@ -526,10 +529,9 class ParametersProc(ProcessingUnit):
526 pairsList.append(pairx)
529 pairsList.append(pairx)
527 pairsList.append(pairy)
530 pairsList.append(pairy)
528
531
529 meteorOps = MeteorOperations()
530 jph = numpy.array([0,0,0,0])
532 jph = numpy.array([0,0,0,0])
531 h = (hmin,hmax)
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 # #Calculate AOA (Error N 3, 4)
536 # #Calculate AOA (Error N 3, 4)
535 # #JONES ET AL. 1998
537 # #JONES ET AL. 1998
@@ -552,10 +554,12 class ParametersProc(ProcessingUnit):
552
554
553 if arrayParameters == None:
555 if arrayParameters == None:
554 self.dataOut.flagNoData = True
556 self.dataOut.flagNoData = True
557 else:
558 self.dataOut.flagNoData = True
555
559
556 return
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 arrayParameters = self.dataOut.data_param
564 arrayParameters = self.dataOut.data_param
561 pairsList = []
565 pairsList = []
@@ -566,12 +570,19 class ParametersProc(ProcessingUnit):
566 jph = numpy.zeros(4)
570 jph = numpy.zeros(4)
567
571
568 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
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 meteorOps = MeteorOperations()
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 h = (hmin,hmax)
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 self.dataOut.data_param = arrayParameters
586 self.dataOut.data_param = arrayParameters
576 return
587 return
577
588
@@ -641,10 +652,11 class ParametersProc(ProcessingUnit):
641 phaseArrival = phaseInt.reshape(phaseInt.size)
652 phaseArrival = phaseInt.reshape(phaseInt.size)
642
653
643 #Dealias
654 #Dealias
644 indAlias = numpy.where(phaseArrival > numpy.pi)
655 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
645 phaseArrival[indAlias] -= 2*numpy.pi
656 # indAlias = numpy.where(phaseArrival > numpy.pi)
646 indAlias = numpy.where(phaseArrival < -numpy.pi)
657 # phaseArrival[indAlias] -= 2*numpy.pi
647 phaseArrival[indAlias] += 2*numpy.pi
658 # indAlias = numpy.where(phaseArrival < -numpy.pi)
659 # phaseArrival[indAlias] += 2*numpy.pi
648
660
649 return phaseDiff, phaseArrival
661 return phaseDiff, phaseArrival
650
662
@@ -2293,17 +2305,25 class PhaseCalibration(Operation):
2293
2305
2294 return False
2306 return False
2295
2307
2296 def __getGammas(self, pairs, k, d, phases):
2308 def __getGammas(self, pairs, d, phases):
2297 gammas = numpy.zeros(2)
2309 gammas = numpy.zeros(2)
2298
2310
2299 for i in range(len(pairs)):
2311 for i in range(len(pairs)):
2300
2312
2301 pairi = pairs[i]
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 #Calculating gamma
2319 #Calculating gamma
2304 jdcos = phases[:,pairi[1]]/(k*d[pairi[1]])
2320 # jdcos = alp1/(k*d1)
2305 jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]])))
2321 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2306
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 #Revised distribution
2327 #Revised distribution
2308 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
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 def __gauss_function(self, t, coeffs):
2366 def __gauss_function(self, t, coeffs):
2347
2367
2348 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2368 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2349
2369
2350 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2370 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2351 meteorOps = MeteorOperations()
2371 meteorOps = MeteorOperations()
2352 nchan = 4
2372 nchan = 4
@@ -2379,14 +2399,14 class PhaseCalibration(Operation):
2379 for iy in range(int(nstepsy)):
2399 for iy in range(int(nstepsy)):
2380 for ix in range(int(nstepsx)):
2400 for ix in range(int(nstepsx)):
2381 jph[pairy[1]] = alpha_y[iy]
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 jph[pairx[1]] = alpha_x[ix]
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 jph_array[:,ix,iy] = jph
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 error = meteorsArray1[:,-1]
2410 error = meteorsArray1[:,-1]
2391 ind1 = numpy.where(error==0)[0]
2411 ind1 = numpy.where(error==0)[0]
2392 penalty[ix,iy] = ind1.size
2412 penalty[ix,iy] = ind1.size
@@ -2402,7 +2422,7 class PhaseCalibration(Operation):
2402 return phOffset
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 dataOut.flagNoData = True
2427 dataOut.flagNoData = True
2408 self.__dataReady = False
2428 self.__dataReady = False
@@ -2435,7 +2455,14 class PhaseCalibration(Operation):
2435 azimuth = 0
2455 azimuth = 0
2436 h = (hmin, hmax)
2456 h = (hmin, hmax)
2437 pairs = ((0,1),(2,3))
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 meteorsArray = self.__buffer
2467 meteorsArray = self.__buffer
2441 error = meteorsArray[:,-1]
2468 error = meteorsArray[:,-1]
@@ -2446,7 +2473,7 class PhaseCalibration(Operation):
2446 phases = meteorsArray[:,8:12]
2473 phases = meteorsArray[:,8:12]
2447
2474
2448 #Calculate Gammas
2475 #Calculate Gammas
2449 gammas = self.__getGammas(pairs, k, distances, phases)
2476 gammas = self.__getGammas(pairs, distances, phases)
2450 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2477 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2451 #Calculate Phases
2478 #Calculate Phases
2452 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
2479 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
@@ -2464,7 +2491,7 class MeteorOperations():
2464
2491
2465 return
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 arrayParameters = arrayParameters0.copy()
2496 arrayParameters = arrayParameters0.copy()
2470 hmin = h[0]
2497 hmin = h[0]
@@ -2476,7 +2503,7 class MeteorOperations():
2476 error = arrayParameters[:,-1]
2503 error = arrayParameters[:,-1]
2477 phases = -arrayParameters[:,8:12] + jph
2504 phases = -arrayParameters[:,8:12] + jph
2478 # phases = numpy.unwrap(phases)
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 #Calculate Heights (Error N 13 and 14)
2508 #Calculate Heights (Error N 13 and 14)
2482 error = arrayParameters[:,-1]
2509 error = arrayParameters[:,-1]
@@ -2491,10 +2518,10 class MeteorOperations():
2491
2518
2492 return arrayParameters
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 arrayAOA = numpy.zeros((phases.shape[0],3))
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 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2526 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2500 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2527 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
@@ -2514,7 +2541,7 class MeteorOperations():
2514 error[indInvalid] = 4
2541 error[indInvalid] = 4
2515 return arrayAOA, error
2542 return arrayAOA, error
2516
2543
2517 def __getDirectionCosines(self, arrayPhase, pairsList):
2544 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2518
2545
2519 #Initializing some variables
2546 #Initializing some variables
2520 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
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 for i in range(2):
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 #First Estimation
2564 #First Estimation
2529 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2565 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
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)
2537
2566
2538 #Most-Accurate Second Estimation
2567 #Most-Accurate Second Estimation
2539 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2568 phi1_aux = ph0 - ph1
2540 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2569 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2541 #Direction Cosine 1
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 #Searching the correct Direction Cosine
2573 #Searching the correct Direction Cosine
2545 cosdir0_aux = cosdir0[:,i]
2574 cosdir0_aux = cosdir0[:,i]
@@ -2557,7 +2586,7 class MeteorOperations():
2557 cosdirY = cosdir[:,1]
2586 cosdirY = cosdir[:,1]
2558
2587
2559 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
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 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2590 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2562
2591
2563 return angles
2592 return angles
@@ -2596,4 +2625,88 class MeteorOperations():
2596 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2625 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2597 error[indInvalid1] = 13
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