schain Commit - r842:ae31287a1940 · Jicamarca Repository

111

self.__updateObjFromInput()

111

self.__updateObjFromInput()

112

self.dataOut.utctimeInit = self.dataIn.utctime

112

self.dataOut.utctimeInit = self.dataIn.utctime

113

self.dataOut.paramInterval = self.dataIn.timeInterval

113

self.dataOut.paramInterval = self.dataIn.timeInterval

114

115

return

114

116

115

#------------------- Get Moments ----------------------------------

117

class SpectralMoments(Operation):

116

def GetMoments(self, channelList = None):

118

117

'''

119

'''

118

Function ~~Get~~Moments()

120

Function SpectralMoments()

121

122

Calculates moments (power, mean, standard deviation) and SNR of the signal

123

124

Type of dataIn: Spectra

119

125

120

Input:

126

Input:

121

channelList : simple channel list to select e.g. [2,3,7]

127

channelList : simple channel list to select e.g. [2,3,7]

122

self.dataOut.data_pre

128

self.dataOut.data_pre : Spectral data

123

self.dataOut.abscissaList

129

self.dataOut.abscissaList : List of frequencies

124

self.dataOut.noise

130

self.dataOut.noise : Noise level per channel

125

131

126

Affected:

132

Affected:

127

self.dataOut.data_param

133

self.dataOut.data_param : Parameters per channel

128

self.dataOut.data_SNR

134

self.dataOut.data_SNR : SNR per channel

129

135

130

'''

136

'''

131

self.dataOut.data_pre = self.dataOut.data_pre[0]

137

132

data = self.dataOut.data_pre

138

def run(self, dataOut, channelList=None):

133

absc = self.dataOut.abscissaList[:-1]

139

134

noise = self.dataOut.noise

140

dataOut.data_pre = dataOut.data_pre[0]

141

data = dataOut.data_pre

142

absc = dataOut.abscissaList[:-1]

143

noise = dataOut.noise

135

144

136

data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))

145

data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))

137

146

138

if channelList== None:

147

if channelList== None:

139

channelList = ~~self~~.~~dataIn~~.channelList

148

channelList = dataOut.channelList

140

~~self~~.dataOut.channelList = channelList

149

dataOut.channelList = channelList

141

150

142

for ind in channelList:

151

for ind in channelList:

143

data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])

152

data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])

144

153

145

~~self~~.dataOut.data_param = data_param[:,1:,:]

154

dataOut.data_param = data_param[:,1:,:]

146

~~self~~.dataOut.data_SNR = data_param[:,0]

155

dataOut.data_SNR = data_param[:,0]

147

return

156

return

148

157

149

def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):

158

def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):

221

#------------------ Get SA Parameters --------------------------

230

#------------------ Get SA Parameters --------------------------

222

231

223

def GetSAParameters(self):

232

def GetSAParameters(self):

233

#SA en frecuencia

224

pairslist = self.dataOut.groupList

234

pairslist = self.dataOut.groupList

225

num_pairs = len(pairslist)

235

num_pairs = len(pairslist)

226

236

249

cspc_par[i,:,:] = __calculateMoments(norm_cspectra)

259

cspc_par[i,:,:] = __calculateMoments(norm_cspectra)

250

#------------------- Get Lags ----------------------------------

260

#------------------- Get Lags ----------------------------------

251

261

252

def GetLags(self):

262

class SALags(Operation):

253

'''

263

'''

254

Function GetMoments()

264

Function GetMoments()

255

265

256

Input:

266

Input:

257

self.dataOut.data_pre

267

self.dataOut.data_pre

258

self.dataOut.abscissaList

268

self.dataOut.abscissaList

259

self.dataOut.noise

269

self.dataOut.noise

260

self.dataOut.normFactor

270

self.dataOut.normFactor

261

self.dataOut.data_SNR

271

self.dataOut.data_SNR

262

self.dataOut.groupList

272

self.dataOut.groupList

263

self.dataOut.nChannels

273

self.dataOut.nChannels

264

274

265

Affected:

275

Affected:

266

self.dataOut.data_param

276

self.dataOut.data_param

267

277

268

'''

278

'''

269

279

def run(self, dataOut):

270

data = ~~self~~.dataOut.data_pre

280

data = dataOut.data_pre

271

normFactor = self.dataOut.normFactor

281

data_acf = dataOut.data_pre[0]

272

nHeights = self.dataOut.nHeights

282

data_ccf = dataOut.data_pre[1]

273

absc = self.dataOut.abscissaList[:-1]

283

274

no~~ise~~ = ~~self~~.dataOut.no~~ise~~

284

normFactor = dataOut.normFactor

275

SNR = self.dataOut.data_SNR

285

nHeights = dataOut.nHeights

276

~~pairsList~~ = ~~self~~.dataOut.~~groupList~~

286

absc = dataOut.abscissaList[:-1]

277

n~~Channels~~ = ~~self~~.dataOut.n~~Channels~~

287

noise = dataOut.noise

278

pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

288

SNR = dataOut.data_SNR

279

self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))

289

# pairsList = dataOut.groupList

290

nChannels = dataOut.nChannels

291

# pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

292

dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))

280

293

281

dataNorm = numpy.abs(data)

294

dataNorm = numpy.abs(data)

282

for l in range(len(pairsList)):

295

for l in range(len(pairsList)):

283

dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]

296

dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]

284

297

285

self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)

298

self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)

286

self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, ~~pairsAutoCorr~~, absc)

299

self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)

287

return

300

return

288

301

289

def __getPairsAutoCorr(self, pairsList, nChannels):

302

# def __getPairsAutoCorr(self, pairsList, nChannels):

290

303

#

291

pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

304

# pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

292

305

#

293

for l in range(len(pairsList)):

306

# for l in range(len(pairsList)):

294

firstChannel = pairsList[l][0]

307

# firstChannel = pairsList[l][0]

295

secondChannel = pairsList[l][1]

308

# secondChannel = pairsList[l][1]

296

309

#

297

#Obteniendo pares de Autocorrelacion

310

# #Obteniendo pares de Autocorrelacion

298

if firstChannel == secondChannel:

311

# if firstChannel == secondChannel:

299

pairsAutoCorr[firstChannel] = int(l)

312

# pairsAutoCorr[firstChannel] = int(l)

300

313

#

301

pairsAutoCorr = pairsAutoCorr.astype(int)

314

# pairsAutoCorr = pairsAutoCorr.astype(int)

302

315

#

303

pairsCrossCorr = range(len(pairsList))

316

# pairsCrossCorr = range(len(pairsList))

304

pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

317

# pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

305

318

#

306

return pairsAutoCorr, pairsCrossCorr

319

# return pairsAutoCorr, pairsCrossCorr

307

320

308

def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):

321

def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):

309

322

331

344

332

return tau

345

return tau

333

346

334

def __calculateLag1Phase(self, data, ~~pairs~~, lagTRange):

347

def __calculateLag1Phase(self, data, lagTRange):

335

data1 = stats.nanmean(data[~~pairs,:,:]~~, axis = 0)

348

data1 = stats.nanmean(data, axis = 0)

336

lag1 = numpy.where(lagTRange == 0)[0][0] + 1

349

lag1 = numpy.where(lagTRange == 0)[0][0] + 1

337

350

338

phase = numpy.angle(data1[lag1,:])

351

phase = numpy.angle(data1[lag1,:])

339

352

340

return phase

353

return phase

341

#------------------- Detect Meteors ------------------------------

342

354

343

def MeteorDetection(self, hei_ref = None, tauindex = 0,

355

class SpectralFitting(Operation):

344

phaseOffsets = None,

356

'''

345

cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,

357

Function GetMoments()

346

noise_timeStep = 4, noise_multiple = 4,

347

multDet_timeLimit = 1, multDet_rangeLimit = 3,

348

phaseThresh = 20, SNRThresh = 5,

349

hmin = 50, hmax=150, azimuth = 0,

350

# channel25X = 0, channel20X = 4, channelCentX = 2,

351

# channel25Y = 1, channel20Y = 3, channelCentY = 2,

352

channelPositions = None) :

353

358

354

'''

355

Function DetectMeteors()

356

Project developed with paper:

357

HOLDSWORTH ET AL. 2004

358

359

Input:

359

Input:

360

self.dataOut.data_pre

360

Output:

361

Variables modified:

362

centerReceiverIndex: From the channels, which is the center receiver

362

'''

363

364

hei_ref: Height reference for the Beacon signal extraction

364

def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):

365

tauindex:

366

predefinedPhaseShifts: Predefined phase offset for the voltge signals

367

368

cohDetection: Whether to user Coherent detection or not

369

cohDet_timeStep: Coherent Detection calculation time step

370

cohDet_thresh: Coherent Detection phase threshold to correct phases

371

372

noise_timeStep: Noise calculation time step

373

noise_multiple: Noise multiple to define signal threshold

374

375

multDet_timeLimit: Multiple Detection Removal time limit in seconds

376

multDet_rangeLimit: Multiple Detection Removal range limit in km

377

378

phaseThresh: Maximum phase difference between receiver to be consider a meteor

379

SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor

380

381

hmin: Minimum Height of the meteor to use it in the further wind estimations

382

hmax: Maximum Height of the meteor to use it in the further wind estimations

383

azimuth: Azimuth angle correction

384

385

Affected:

386

self.dataOut.data_param

387

388

Rejection Criteria (Errors):

389

0: No error; analysis OK

390

1: SNR < SNR threshold

391

2: angle of arrival (AOA) ambiguously determined

392

3: AOA estimate not feasible

393

4: Large difference in AOAs obtained from different antenna baselines

394

5: echo at start or end of time series

395

6: echo less than 5 examples long; too short for analysis

396

7: echo rise exceeds 0.3s

397

8: echo decay time less than twice rise time

398

9: large power level before echo

399

10: large power level after echo

400

11: poor fit to amplitude for estimation of decay time

401

12: poor fit to CCF phase variation for estimation of radial drift velocity

402

13: height unresolvable echo: not valid height within 70 to 110 km

403

14: height ambiguous echo: more then one possible height within 70 to 110 km

404

15: radial drift velocity or projected horizontal velocity exceeds 200 m/s

405

16: oscilatory echo, indicating event most likely not an underdense echo

406

407

17: phase difference in meteor Reestimation

408

409

Data Storage:

410

Meteors for Wind Estimation (8):

411

Utc Time | Range Height

412

Azimuth Zenith errorCosDir

413

VelRad errorVelRad

414

Phase0 Phase1 Phase2 Phase3

415

TypeError

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

365

425

#Get Beacon signal

426

newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])

427

366

428

if h~~ei_ref~~ != None:

367

if path != None:

429

newheis = numpy.where(self.dataOut.heightList>hei_ref)

368

sys.path.append(path)

369

self.dataOut.library = importlib.import_module(file)

430

370

431

heiRang = self.dataOut.getHeiRange()

371

#To be inserted as a parameter

432

372

groupArray = numpy.array(groupList)

433

# nChannel = self.dataOut.nChannels

373

# groupArray = numpy.array([[0,1],[2,3]])

434

# for i in range(nChannel):

374

self.dataOut.groupList = groupArray

435

# if i != centerReceiverIndex:

436

# pairslist.append((centerReceiverIndex,i))

437

438

#****************REMOVING HARDWARE PHASE DIFFERENCES***************

439

# see if the user put in pre defined phase shifts

440

voltsPShift = self.dataOut.data_pre.copy()

441

375

442

# if predefinedPhaseShifts != None:

376

nGroups = groupArray.shape[0]

443

# hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180

377

nChannels = self.dataIn.nChannels

444

#

378

nHeights=self.dataIn.heightList.size

445

# # elif beaconPhaseShifts:

446

# # #get hardware phase shifts using beacon signal

447

# # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)

448

# # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)

449

#

450

# else:

451

# hardwarePhaseShifts = numpy.zeros(5)

452

#

453

# voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')

454

# for i in range(self.dataOut.data_pre.shape[0]):

455

# voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])

456

457

#******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********

458

459

#Remove DC

460

voltsDC = numpy.mean(voltsPShift,1)

461

voltsDC = numpy.mean(voltsDC,1)

462

for i in range(voltsDC.shape[0]):

463

voltsPShift[i] = voltsPShift[i] - voltsDC[i]

464

465

#Don't considerate last heights, theyre used to calculate Hardware Phase Shift

466

voltsPShift = voltsPShift[:,:,:newheis[0][0]]

467

379

468

#************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********

380

#Parameters Array

469

#Coherent Detection

381

self.dataOut.data_param = None

470

if cohDetection:

471

#use coherent detection to get the net power

472

cohDet_thresh = cohDet_thresh*numpy.pi/180

473

voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh)

474

382

475

#Non-coherent detection!

383

#Set constants

476

powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)

384

constants = self.dataOut.library.setConstants(self.dataIn)

477

#********** END OF COH/NON-COH POWER CALCULATION**********************

385

self.dataOut.constants = constants

478

386

M = self.dataIn.normFactor

479

#********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************

387

N = self.dataIn.nFFTPoints

480

#Get noise

388

ippSeconds = self.dataIn.ippSeconds

481

noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)

389

K = self.dataIn.nIncohInt

482

# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)

390

pairsArray = numpy.array(self.dataIn.pairsList)

483

#Get signal threshold

484

signalThresh = noise_multiple*noise

485

#Meteor echoes detection

486

listMeteors = self.__findMeteors(powerNet, signalThresh)

487

#******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********

488

391

489

#************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************

392

#List of possible combinations

490

#Parameters

393

listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)

491

heiRange = self.dataOut.getHeiRange()

394

indCross = numpy.zeros(len(list(listComb)), dtype = 'int')

492

rangeInterval = heiRange[1] - heiRange[0]

493

rangeLimit = multDet_rangeLimit/rangeInterval

494

timeLimit = multDet_timeLimit/self.dataOut.timeInterval

495

#Multiple detection removals

496

listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)

497

#************ END OF REMOVE MULTIPLE DETECTIONS **********************

498

395

499

#********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************

396

if getSNR:

500

#Parameters

397

listChannels = groupArray.reshape((groupArray.size))

501

phaseThresh = phaseThresh*numpy.pi/180

398

listChannels.sort()

502

thresh = [phaseThresh, noise_multiple, SNRThresh]

399

noise = self.dataIn.getNoise()

503

#Meteor reestimation (Errors N 1, 6, 12, 17)

400

self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])

504

listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)

505

# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)

506

#Estimation of decay times (Errors N 7, 8, 11)

507

listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)

508

#******************* END OF METEOR REESTIMATION *******************

509

401

510

#********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************

402

for i in range(nGroups):

511

#Calculating Radial Velocity (Error N 15)

403

coord = groupArray[i,:]

512

radialStdThresh = 10

513

listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval)

514

515

if len(listMeteors4) > 0:

516

#Setting New Array

517

date = self.dataOut.utctime

518

arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)

519

520

#Correcting phase offset

521

if phaseOffsets != None:

522

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

523

arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

524

404

525

#~~Second Pairslist~~

405

#Input data array

526

pairsList = []

406

data = self.dataIn.data_spc[coord,:,:]/(M*N)

527

pairx = (0,1)

407

data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))

528

pairy = (2,3)

529

pairsList.append(pairx)

530

pairsList.append(pairy)

531

408

532

jph = numpy.array([0,0,0,0])

409

#Cross Spectra data array for Covariance Matrixes

533

h = (~~hmin~~,~~hmax~~)

410

ind = 0

534

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

411

for pairs in listComb:

412

pairsSel = numpy.array([coord[x],coord[y]])

413

indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])

414

ind += 1

415

dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)

416

dataCross = dataCross**2/K

535

417

536

# #Calculate AOA (Error N 3, 4)

418

for h in range(nHeights):

537

# #JONES ET AL. 1998

419

# print self.dataOut.heightList[h]

538

# error = arrayParameters[:,-1]

539

# AOAthresh = numpy.pi/8

540

# phases = -arrayParameters[:,9:13]

541

# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)

542

#

543

# #Calculate Heights (Error N 13 and 14)

544

# error = arrayParameters[:,-1]

545

# Ranges = arrayParameters[:,2]

546

# zenith = arrayParameters[:,5]

547

# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)

548

# error = arrayParameters[:,-1]

549

#********************* END OF PARAMETERS CALCULATION **************************

550

420

551

#***************************+ PASS DATA TO NEXT STEP **********************

421

#Input

552

# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))

422

d = data[:,h]

553

self.dataOut.data_param = arrayParameters

423

554

424

#Covariance Matrix

555

if arrayParameters == None:

425

D = numpy.diag(d**2/K)

556

self.dataOut.flagNoData = True

426

ind = 0

557

else:

427

for pairs in listComb:

558

self.dataOut.flagNoData = True

428

#Coordinates in Covariance Matrix

559

429

x = pairs[0]

430

y = pairs[1]

431

#Channel Index

432

S12 = dataCross[ind,:,h]

433

D12 = numpy.diag(S12)

434

#Completing Covariance Matrix with Cross Spectras

435

D[x*N:(x+1)*N,y*N:(y+1)*N] = D12

436

D[y*N:(y+1)*N,x*N:(x+1)*N] = D12

437

ind += 1

438

Dinv=numpy.linalg.inv(D)

439

L=numpy.linalg.cholesky(Dinv)

440

LT=L.T

441

442

dp = numpy.dot(LT,d)

443

444

#Initial values

445

data_spc = self.dataIn.data_spc[coord,:,h]

446

447

if (h>0)and(error1[3]<5):

448

p0 = self.dataOut.data_param[i,:,h-1]

449

else:

450

p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))

451

452

try:

453

#Least Squares

454

minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)

455

# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))

456

#Chi square error

457

error0 = numpy.sum(infodict['fvec']**2)/(2*N)

458

#Error with Jacobian

459

error1 = self.dataOut.library.errorFunction(minp,constants,LT)

460

except:

461

minp = p0*numpy.nan

462

error0 = numpy.nan

463

error1 = p0*numpy.nan

464

465

#Save

466

if self.dataOut.data_param == None:

467

self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan

468

self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan

469

470

self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))

471

self.dataOut.data_param[i,:,h] = minp

560

return

472

return

561

473

562

def CorrectMeteorPhases(self, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):

474

def __residFunction(self, p, dp, LT, constants):

563

564

arrayParameters = self.dataOut.data_param

565

pairsList = []

566

pairx = (0,1)

567

pairy = (2,3)

568

pairsList.append(pairx)

569

pairsList.append(pairy)

570

jph = numpy.zeros(4)

571

572

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

573

# arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

574

arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))

575

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)

582

h = (hmin,hmax)

583

475

584

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

476

fm = self.dataOut.library.modelFunction(p, constants)

477

fmp=numpy.dot(LT,fm)

478

479

return dp-fmp

480

481

def __getSNR(self, z, noise):

585

482

586

self.dataOut.data_param = arrayParameters

483

avg = numpy.average(z, axis=1)

587

return

484

SNR = (avg.T-noise)/noise

485

SNR = SNR.T

486

return SNR

588

487

589

def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):

488

def __chisq(p,chindex,hindex):

590

489

#similar to Resid but calculates CHI**2

591

minIndex = min(newheis[0])

490

[LT,d,fm]=setupLTdfm(p,chindex,hindex)

592

maxIndex = max(newheis[0])

491

dp=numpy.dot(LT,d)

593

492

fmp=numpy.dot(LT,fm)

594

voltage = voltage0[:,:,minIndex:maxIndex+1]

493

chisq=numpy.dot((dp-fmp).T,(dp-fmp))

595

nLength = voltage.shape[1]/n

494

return chisq

596

nMin = 0

495

597

nMax = 0

496

class WindProfiler(Operation):

598

phaseOffset = numpy.zeros((len(pairslist),n))

497

599

498

__isConfig = False

600

for i in range(n):

601

nMax += nLength

602

phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))

603

phaseCCF = numpy.mean(phaseCCF, axis = 2)

604

phaseOffset[:,i] = phaseCCF.transpose()

605

nMin = nMax

606

# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)

607

608

#Remove Outliers

609

factor = 2

610

wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))

611

dw = numpy.std(wt,axis = 1)

612

dw = dw.reshape((dw.size,1))

613

ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))

614

phaseOffset[ind] = numpy.nan

615

phaseOffset = stats.nanmean(phaseOffset, axis=1)

616

499

617

return phaseOffset

500

__initime = None

501

__lastdatatime = None

502

__integrationtime = None

618

503

619

def __shiftPhase(self, data, phaseShift):

504

__buffer = None

620

#this will shift the phase of a complex number

621

dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)

622

return dataShifted

623

505

624

def __estimatePhaseDifference(self, array, pairslist):

506

__dataReady = False

625

nChannel = array.shape[0]

507

626

nHeights = array.shape[2]

508

__firstdata = None

627

numPairs = len(pairslist)

509

628

# phaseCCF = numpy.zeros((nChannel, 5, nHeights))

510

n = None

629

phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))

511

512

def __init__(self):

513

Operation.__init__(self)

514

515

def __calculateCosDir(self, elev, azim):

516

zen = (90 - elev)*numpy.pi/180

517

azim = azim*numpy.pi/180

518

cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))

519

cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)

630

520

631

#Correct phases

521

signX = numpy.sign(numpy.cos(azim))

632

derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]

522

signY = numpy.sign(numpy.sin(azim))

633

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

634

523

635

if indDer[0].shape[0] > 0:

524

cosDirX = numpy.copysign(cosDirX, signX)

636

for i in range(indDer[0].shape[0]):

525

cosDirY = numpy.copysign(cosDirY, signY)

637

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])

526

return cosDirX, cosDirY

638

phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi

639

527

640

# for j in range(numSides):

528

def __calculateAngles(self, theta_x, theta_y, azimuth):

641

# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])

529

642

# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)

530

dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)

643

#

531

zenith_arr = numpy.arccos(dir_cosw)

644

#Linear

532

azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180

645

phaseInt = numpy.zeros((numPairs,1))

646

angAllCCF = phaseCCF[:,[0,1,3,4],0]

647

for j in range(numPairs):

648

fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])

649

phaseInt[j] = fit[1]

650

#Phase Differences

651

phaseDiff = phaseInt - phaseCCF[:,2,:]

652

phaseArrival = phaseInt.reshape(phaseInt.size)

653

533

654

#Dealias

534

dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)

655

phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))

535

dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)

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

660

536

661

return phaseDiff, phaseArrival

537

return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw

662

538

663

def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):

539

def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):

664

#this function will run the coherent detection used in Holdworth et al. 2004 and return the net power

665

#find the phase shifts of each channel over 1 second intervals

666

#only look at ranges below the beacon signal

667

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

668

numBlocks = int(volts.shape[1]/numProfPerBlock)

669

numHeights = volts.shape[2]

670

nChannel = volts.shape[0]

671

voltsCohDet = volts.copy()

672

540

673

pairsarray = numpy.array(pairslist)

674

indSides = pairsarray[:,1]

675

# indSides = numpy.array(range(nChannel))

676

# indSides = numpy.delete(indSides, indCenter)

677

#

541

#

678

# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)

542

if horOnly:

679

listBlocks = numpy.array_split(volts, numBlocks, 1)

543

A = numpy.c_[dir_cosu,dir_cosv]

544

else:

545

A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]

546

A = numpy.asmatrix(A)

547

A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()

548

549

return A1

550

551

def __correctValues(self, heiRang, phi, velRadial, SNR):

552

listPhi = phi.tolist()

553

maxid = listPhi.index(max(listPhi))

554

minid = listPhi.index(min(listPhi))

680

555

681

startInd = 0

556

rango = range(len(phi))

682

endInd = 0

557

# rango = numpy.delete(rango,maxid)

558

559

heiRang1 = heiRang*math.cos(phi[maxid])

560

heiRangAux = heiRang*math.cos(phi[minid])

561

indOut = (heiRang1 < heiRangAux[0]).nonzero()

562

heiRang1 = numpy.delete(heiRang1,indOut)

563

564

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

565

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

566

567

for i in rango:

568

x = heiRang*math.cos(phi[i])

569

y1 = velRadial[i,:]

570

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

683

571

684

for i in range(numBlocks):

572

x1 = heiRang1

685

startInd = endInd

573

y11 = f1(x1)

686

endInd = endInd + listBlocks[i].shape[1]

687

574

688

~~arrayBlock~~ = ~~listBlocks~~[i]

575

y2 = SNR[i,:]

689

# arrayBlockCenter = listCenter[i]

576

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

577

y21 = f2(x1)

690

578

691

#Estimate the Phase Difference

579

velRadial1[i,:] = y11

692

phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)

580

SNR1[i,:] = y21

693

#Phase Difference RMS

581

694

arrayPhaseRMS = numpy.abs(phaseDiff)

582

return heiRang1, velRadial1, SNR1

695

phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)

583

696

indPhase = numpy.where(phaseRMSaux==4)

584

def __calculateVelUVW(self, A, velRadial):

697

#Shifting

698

if indPhase[0].shape[0] > 0:

699

for j in range(indSides.size):

700

arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())

701

voltsCohDet[:,startInd:endInd,:] = arrayBlock

702

703

return voltsCohDet

704

705

def __calculateCCF(self, volts, pairslist ,laglist):

706

585

707

nHeights = volts.shape[2]

586

#Operacion Matricial

708

nPoints = volts.shape[1]

587

# velUVW = numpy.zeros((velRadial.shape[1],3))

709

voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')

588

# for ind in range(velRadial.shape[1]):

589

# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])

590

# velUVW = velUVW.transpose()

591

velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))

592

velUVW[:,:] = numpy.dot(A,velRadial)

710

593

711

for i in range(len(pairslist)):

594

712

volts1 = volts[pairslist[i][0]]

595

return velUVW

713

volts2 = volts[pairslist[i][1]]

714

715

for t in range(len(laglist)):

716

idxT = laglist[t]

717

if idxT >= 0:

718

vStacked = numpy.vstack((volts2[idxT:,:],

719

numpy.zeros((idxT, nHeights),dtype='complex')))

720

else:

721

vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),

722

volts2[:(nPoints + idxT),:]))

723

voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)

724

596

725

vStacked = None

597

# def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):

726

return voltsCCF

598

def techniqueDBS(self, velRadial0, heiRang, SNR0, kwargs):

599

"""

600

Function that implements Doppler Beam Swinging (DBS) technique.

727

601

728

def __getNoise(self, power, timeSegment, timeInterval):

602

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

729

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

603

Direction correction (if necessary), Ranges and SNR

730

numBlocks = int(power.shape[0]/numProfPerBlock)

731

numHeights = power.shape[1]

732

733

listPower = numpy.array_split(power, numBlocks, 0)

734

noise = numpy.zeros((power.shape[0], power.shape[1]))

735

noise1 = numpy.zeros((power.shape[0], power.shape[1]))

736

604

737

startInd = 0

605

Output: Winds estimation (Zonal, Meridional and Vertical)

738

endInd = 0

739

606

740

for i in range(numBlocks): #split por canal

607

Parameters affected: Winds, height range, SNR

741

startInd = endInd

608

"""

742

endInd = endInd + listPower[i].shape[0]

743

744

arrayBlock = listPower[i]

745

noiseAux = numpy.mean(arrayBlock, 0)

746

# noiseAux = numpy.median(noiseAux)

747

# noiseAux = numpy.mean(arrayBlock)

748

noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux

749

750

noiseAux1 = numpy.mean(arrayBlock)

751

noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1

752

753

return noise, noise1

754

755

def __findMeteors(self, power, thresh):

756

nProf = power.shape[0]

757

nHeights = power.shape[1]

758

listMeteors = []

759

609

760

for i in range(nHeights):

610

if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):

761

powerAux = power[:,i]

611

theta_x = numpy.array(kwargs['dirCosx'])

762

threshAux = thresh[:,i]

612

theta_y = numpy.array(kwargs['dirCosy'])

763

613

else:

764

indUPthresh = numpy.where(powerAux > threshAux)[0]

614

elev = numpy.array(kwargs['elevation'])

765

indDNthresh = numpy.where(powerAux <= threshAux)[0]

615

azim = numpy.array(kwargs['azimuth'])

766

616

theta_x, theta_y = self.__calculateCosDir(elev, azim)

767

j = 0

617

azimuth = kwargs['correctAzimuth']

768

618

if kwargs.has_key('horizontalOnly'):

769

while (j < indUPthresh.size - 2):

619

horizontalOnly = kwargs['horizontalOnly']

770

if (indUPthresh[j + 2] == indUPthresh[j] + 2):

620

else: horizontalOnly = False

771

indDNAux = numpy.where(indDNthresh > indUPthresh[j])

621

if kwargs.has_key('correctFactor'):

772

indDNthresh = indDNthresh[indDNAux]

622

correctFactor = kwargs['correctFactor']

773

623

else: correctFactor = 1

774

if (indDNthresh.size > 0):

624

if kwargs.has_key('channelList'):

775

indEnd = indDNthresh[0] - 1

625

channelList = kwargs['channelList']

776

indInit = indUPthresh[j]

626

if len(channelList) == 2:

777

627

horizontalOnly = True

778

meteor = powerAux[indInit:indEnd + 1]

628

arrayChannel = numpy.array(channelList)

779

indPeak = meteor.argmax() + indInit

629

param = param[arrayChannel,:,:]

780

FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))

630

theta_x = theta_x[arrayChannel]

781

631

theta_y = theta_y[arrayChannel]

782

listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!

632

783

j = numpy.where(indUPthresh == indEnd)[0] + 1

633

azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)

784

else: j+=1

634

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)

785

else: j+=1

635

A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)

786

636

787

return listMeteors

637

#Calculo de Componentes de la velocidad con DBS

638

winds = self.__calculateVelUVW(A,velRadial1)

639

640

return winds, heiRang1, SNR1

788

641

789

def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):

642

def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):

790

643

791

~~arrayMeteors~~ = numpy.asarray(~~listMeteors~~)

644

posx = numpy.asarray(posx)

792

listMeteors1 = []

645

posy = numpy.asarray(posy)

793

646

794

while arrayMeteors.shape[0] > 0:

647

#Rotacion Inversa para alinear con el azimuth

795

FLAs = arrayMeteors[:,4]

648

if azimuth!= None:

796

maxFLA = FLAs.argmax()

649

azimuth = azimuth*math.pi/180

797

listMeteors1.append(arrayMeteors[maxFLA,:])

650

posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)

798

651

posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)

799

MeteorInitTime = arrayMeteors[maxFLA,1]

652

else:

800

MeteorEndTime = arrayMeteors[maxFLA,3]

653

posx1 = posx

801

MeteorHeight = arrayMeteors[maxFLA,0]

654

posy1 = posy

802

655

803

#Check neighborhood

656

#Calculo de Distancias

804

maxHeightIndex = MeteorHeight + rangeLimit

657

distx = numpy.zeros(pairsCrossCorr.size)

805

minHeightIndex = MeteorHeight - rangeLimit

658

disty = numpy.zeros(pairsCrossCorr.size)

806

minTimeIndex = MeteorInitTime - timeLimit

659

dist = numpy.zeros(pairsCrossCorr.size)

807

maxTimeIndex = MeteorEndTime + timeLimit

660

ang = numpy.zeros(pairsCrossCorr.size)

808

661

809

#Check Heights

662

for i in range(pairsCrossCorr.size):

810

indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)

663

distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]

811

indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)

664

disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]

812

indBoth = numpy.where(numpy.logical_and(indTime,indHeight))

665

dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)

666

ang[i] = numpy.arctan2(disty[i],distx[i])

667

#Calculo de Matrices

668

nPairs = len(pairs)

669

ang1 = numpy.zeros((nPairs, 2, 1))

670

dist1 = numpy.zeros((nPairs, 2, 1))

671

672

for j in range(nPairs):

673

dist1[j,0,0] = dist[pairs[j][0]]

674

dist1[j,1,0] = dist[pairs[j][1]]

675

ang1[j,0,0] = ang[pairs[j][0]]

676

ang1[j,1,0] = ang[pairs[j][1]]

813

677

814

arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)

678

return distx,disty, dist1,ang1

679

680

def __calculateVelVer(self, phase, lagTRange, _lambda):

681

682

Ts = lagTRange[1] - lagTRange[0]

683

velW = -_lambda*phase/(4*math.pi*Ts)

815

684

816

return ~~listMeteors1~~

685

return velW

817

686

818

def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):

687

def __calculateVelHorDir(self, dist, tau1, tau2, ang):

819

n~~umHeight~~s = ~~volts~~.shape[2]

688

nPairs = tau1.shape[0]

820

nChannel = volts.shape[0]

689

vel = numpy.zeros((nPairs,3,tau1.shape[2]))

821

690

822

thresholdPhase = thresh[0]

691

angCos = numpy.cos(ang)

823

thresholdNoise = thresh[1]

692

angSin = numpy.sin(ang)

824

thresholdDB = float(thresh[2])

825

693

826

thresholdDB1 = 10**(thresholdDB/10)

694

vel0 = dist*tau1/(2*tau2**2)

827

pairsarray = numpy.array(pairslist)

695

vel[:,0,:] = (vel0*angCos).sum(axis = 1)

828

indSides = pairsarray[:,1]

696

vel[:,1,:] = (vel0*angSin).sum(axis = 1)

829

697

830

pairslist1 = list(pairslist)

698

ind = numpy.where(numpy.isinf(vel))

831

pairslist1.append((0,1))

699

vel[ind] = numpy.nan

832

pairslist1.append((3,4))

700

701

return vel

702

703

def __getPairsAutoCorr(self, pairsList, nChannels):

833

704

834

listMeteors1 = []

705

pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

835

listPowerSeries = []

706

836

listVoltageSeries = []

707

for l in range(len(pairsList)):

837

#volts has the war data

708

firstChannel = pairsList[l][0]

709

secondChannel = pairsList[l][1]

710

711

#Obteniendo pares de Autocorrelacion

712

if firstChannel == secondChannel:

713

pairsAutoCorr[firstChannel] = int(l)

714

715

pairsAutoCorr = pairsAutoCorr.astype(int)

838

716

839

if frequency == 30e6:

717

pairsCrossCorr = range(len(pairsList))

840

timeLag = 45*10**-3

718

pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

841

else:

842

timeLag = 15*10**-3

843

lag = numpy.ceil(timeLag/timeInterval)

844

719

845

for i in range(len(listMeteors)):

720

return pairsAutoCorr, pairsCrossCorr

846

721

847

###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################

722

def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):

848

meteorAux = numpy.zeros(16)

723

"""

849

724

Function that implements Spaced Antenna (SA) technique.

850

#Loading meteor Data (mHeight, mStart, mPeak, mEnd)

851

mHeight = listMeteors[i][0]

852

mStart = listMeteors[i][1]

853

mPeak = listMeteors[i][2]

854

mEnd = listMeteors[i][3]

855

856

#get the volt data between the start and end times of the meteor

857

meteorVolts = volts[:,mStart:mEnd+1,mHeight]

858

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

859

860

#3.6. Phase Difference estimation

861

phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)

862

863

#3.7. Phase difference removal & meteor start, peak and end times reestimated

864

#meteorVolts0.- all Channels, all Profiles

865

meteorVolts0 = volts[:,:,mHeight]

866

meteorThresh = noise[:,mHeight]*thresholdNoise

867

meteorNoise = noise[:,mHeight]

868

meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting

869

powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power

870

871

#Times reestimation

872

mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]

873

if mStart1.size > 0:

874

mStart1 = mStart1[-1] + 1

875

876

else:

877

mStart1 = mPeak

878

879

mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1

880

mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]

881

if mEndDecayTime1.size == 0:

882

mEndDecayTime1 = powerNet0.size

883

else:

884

mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1

885

# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()

886

887

#meteorVolts1.- all Channels, from start to end

888

meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]

889

meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]

890

if meteorVolts2.shape[1] == 0:

891

meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]

892

meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)

893

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)

894

##################### END PARAMETERS REESTIMATION #########################

895

896

##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################

897

# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis

898

if meteorVolts2.shape[1] > 0:

899

#Phase Difference re-estimation

900

phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation

901

# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)

902

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])

903

phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))

904

meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting

905

906

#Phase Difference RMS

907

phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))

908

powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)

909

#Data from Meteor

910

mPeak1 = powerNet1.argmax() + mStart1

911

mPeakPower1 = powerNet1.max()

912

noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])

913

mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux

914

Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])

915

Meteor1 = numpy.hstack((Meteor1,phaseDiffint))

916

PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]

917

#Vectorize

918

meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]

919

meteorAux[7:11] = phaseDiffint[0:4]

920

921

#Rejection Criterions

922

if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation

923

meteorAux[-1] = 17

924

elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB

925

meteorAux[-1] = 1

926

927

928

else:

929

meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]

930

meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis

931

PowerSeries = 0

932

933

listMeteors1.append(meteorAux)

934

listPowerSeries.append(PowerSeries)

935

listVoltageSeries.append(meteorVolts1)

936

937

return listMeteors1, listPowerSeries, listVoltageSeries

938

939

def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):

940

725

941

threshError = 10

726

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

942

#Depending if it is 30 or 50 MHz

727

Direction correction (if necessary), Ranges and SNR

943

if frequency == 30e6:

944

timeLag = 45*10**-3

945

else:

946

timeLag = 15*10**-3

947

lag = numpy.ceil(timeLag/timeInterval)

948

728

949

listMeteors1 = []

729

Output: Winds estimation (Zonal, Meridional and Vertical)

950

730

951

for i in range(len(listMeteors)):

731

Parameters affected: Winds

952

meteorPower = listPower[i]

732

"""

953

meteorAux = listMeteors[i]

733

#Cross Correlation pairs obtained

954

734

pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

955

if meteorAux[-1] == 0:

735

pairsArray = numpy.array(pairsList)[pairsCrossCorr]

956

736

pairsSelArray = numpy.array(pairsSelected)

957

try:

737

pairs = []

958

indmax = meteorPower.argmax()

959

indlag = indmax + lag

960

961

y = meteorPower[indlag:]

962

x = numpy.arange(0, y.size)*timeLag

963

964

#first guess

965

a = y[0]

966

tau = timeLag

967

#exponential fit

968

popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])

969

y1 = self.__exponential_function(x, *popt)

970

#error estimation

971

error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))

972

973

decayTime = popt[1]

974

riseTime = indmax*timeInterval

975

meteorAux[11:13] = [decayTime, error]

976

977

#Table items 7, 8 and 11

978

if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s

979

meteorAux[-1] = 7

980

elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time

981

meteorAux[-1] = 8

982

if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time

983

meteorAux[-1] = 11

984

985

986

except:

987

meteorAux[-1] = 11

988

989

990

listMeteors1.append(meteorAux)

991

738

992

return listMeteors1

739

#Wind estimation pairs obtained

993

740

for i in range(pairsSelArray.shape[0]/2):

994

#Exponential Function

741

ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]

742

ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]

743

pairs.append((ind1,ind2))

744

745

indtau = tau.shape[0]/2

746

tau1 = tau[:indtau,:]

747

tau2 = tau[indtau:-1,:]

748

tau1 = tau1[pairs,:]

749

tau2 = tau2[pairs,:]

750

phase1 = tau[-1,:]

751

752

#---------------------------------------------------------------------

753

#Metodo Directo

754

distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)

755

winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)

756

winds = stats.nanmean(winds, axis=0)

757

#---------------------------------------------------------------------

758

#Metodo General

759

# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)

760

# #Calculo Coeficientes de Funcion de Correlacion

761

# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)

762

# #Calculo de Velocidades

763

# winds = self.calculateVelUV(F,G,A,B,H)

995

764

996

def __exponential_function(self, x, a, tau):

765

#---------------------------------------------------------------------

997

y = a*numpy.exp(-x/tau)

766

winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)

998

return y

767

winds = correctFactor*winds

768

return winds

999

769

1000

def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):

770

def __checkTime(self, currentTime, paramInterval, outputInterval):

1001

771

1002

pairslist1 = list(pairslist)

772

dataTime = currentTime + paramInterval

1003

pairslist1.append((0,1))

773

deltaTime = dataTime - self.__initime

1004

pairslist1.append((3,4))

1005

numPairs = len(pairslist1)

1006

#Time Lag

1007

timeLag = 45*10**-3

1008

c = 3e8

1009

lag = numpy.ceil(timeLag/timeInterval)

1010

freq = 30e6

1011

774

1012

listMeteors1 = []

775

if deltaTime >= outputInterval or deltaTime < 0:

776

self.__dataReady = True

777

return

778

779

def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):

780

'''

781

Function that implements winds estimation technique with detected meteors.

1013

782

1014

for i in range(len(listMeteors)):

783

Input: Detected meteors, Minimum meteor quantity to wind estimation

1015

meteorAux = listMeteors[i]

784

1016

if meteorAux[-1] == 0:

785

Output: Winds estimation (Zonal and Meridional)

1017

mStart = listMeteors[i][1]

786

1018

mPeak = listMeteors[i][2]

787

Parameters affected: Winds

1019

mLag = mPeak - mStart + lag

788

'''

1020

789

# print arrayMeteor.shape

1021

#get the volt data between the start and end times of the meteor

790

#Settings

1022

meteorVolts = listVolts[i]

791

nInt = (heightMax - heightMin)/2

1023

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

792

# print nInt

1024

793

nInt = int(nInt)

1025

#Get CCF

794

# print nInt

1026

allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])

795

winds = numpy.zeros((2,nInt))*numpy.nan

1027

796

1028

#Method 2

797

#Filter errors

1029

slopes = numpy.zeros(numPairs)

798

error = numpy.where(arrayMeteor[:,-1] == 0)[0]

1030

time = numpy.array([-2,-1,1,2])*timeInterval

799

finalMeteor = arrayMeteor[error,:]

1031

angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])

800

1032

801

#Meteor Histogram

1033

#Correct phases

802

finalHeights = finalMeteor[:,2]

1034

derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]

803

hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))

1035

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

804

nMeteorsPerI = hist[0]

805

heightPerI = hist[1]

806

807

#Sort of meteors

808

indSort = finalHeights.argsort()

809

finalMeteor2 = finalMeteor[indSort,:]

810

811

# Calculating winds

812

ind1 = 0

813

ind2 = 0

814

815

for i in range(nInt):

816

nMet = nMeteorsPerI[i]

817

ind1 = ind2

818

ind2 = ind1 + nMet

819

820

meteorAux = finalMeteor2[ind1:ind2,:]

821

822

if meteorAux.shape[0] >= meteorThresh:

823

vel = meteorAux[:, 6]

824

zen = meteorAux[:, 4]*numpy.pi/180

825

azim = meteorAux[:, 3]*numpy.pi/180

1036

826

1037

if indDer[0].shape[0] > 0:

827

n = numpy.cos(zen)

1038

for i in range(indDer[0].shape[0]):

828

# m = (1 - n**2)/(1 - numpy.tan(azim)**2)

1039

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])

829

# l = m*numpy.tan(azim)

1040

angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi

830

l = numpy.sin(zen)*numpy.sin(azim)

1041

831

m = numpy.sin(zen)*numpy.cos(azim)

1042

# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))

1043

for j in range(numPairs):

1044

fit = stats.linregress(time, angAllCCF[j,:])

1045

slopes[j] = fit[0]

1046

832

1047

#Remove Outlier

833

A = numpy.vstack((l, m)).transpose()

1048

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

834

A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())

1049

# ~~slopes = numpy.delete(slopes,indOut~~)

835

windsAux = numpy.dot(A1, vel)

1050

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

1051

# slopes = numpy.delete(slopes,indOut)

1052

1053

radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)

1054

radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)

1055

meteorAux[-2] = radialError

1056

meteorAux[-3] = radialVelocity

1057

836

1058

#Setting Error

837

winds[0,i] = windsAux[0]

1059

#Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s

838

winds[1,i] = windsAux[1]

1060

if numpy.abs(radialVelocity) > 200:

1061

meteorAux[-1] = 15

1062

#Number 12: Poor fit to CCF variation for estimation of radial drift velocity

1063

elif radialError > radialStdThresh:

1064

meteorAux[-1] = 12

1065

839

1066

listMeteors1.append(meteorAux)

840

return winds, heightPerI[:-1]

1067

return listMeteors1

1068

841

1069

def __setNewArrays(self, listMeteors, date, heiRang):

842

def techniqueNSM_SA(self, **kwargs):

1070

843

metArray = kwargs['metArray']

1071

#New arrays

844

heightList = kwargs['heightList']

1072

arrayMeteors = numpy.array(listMeteors)

845

timeList = kwargs['timeList']

1073

arrayParameters = numpy.zeros((len(listMeteors), 13))

1074

846

1075

#Date inclusion

847

rx_location = kwargs['rx_location']

1076

# date = re.findall(r'\((.*?)\)', date)

848

groupList = kwargs['groupList']

1077

# date = date[0].split(',')

849

azimuth = kwargs['azimuth']

1078

# date = map(int, date)

850

dfactor = kwargs['dfactor']

1079

#

851

k = kwargs['k']

1080

# if len(date)<6:

1081

# date.append(0)

1082

#

1083

# date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]

1084

# arrayDate = numpy.tile(date, (len(listMeteors), 1))

1085

arrayDate = numpy.tile(date, (len(listMeteors)))

1086

852

1087

#Meteor array

853

azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)

1088

# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]

854

d = dist*dfactor

1089

# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))

855

#Phase calculation

856

metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)

1090

857

1091

#Parameters Array

858

metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities

1092

arrayParameters[:,0] = arrayDate #Date

1093

arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range

1094

arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error

1095

arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases

1096

arrayParameters[:,-1] = arrayMeteors[:,-1] #Error

1097

1098

859

1099

return arrayParameters

860

velEst = numpy.zeros((heightList.size,2))*numpy.nan

1100

861

azimuth1 = azimuth1*numpy.pi/180

1101

# def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):

862

1102

#

863

for i in range(heightList.size):

1103

# arrayAOA = numpy.zeros((phases.shape[0],3))

864

h = heightList[i]

1104

# cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)

865

indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]

1105

#

866

metHeight = metArray1[indH,:]

1106

# arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)

867

if metHeight.shape[0] >= 2:

1107

# cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)

868

velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities

1108

# arrayAOA[:,2] = cosDirError

869

iazim = metHeight[:,1].astype(int)

1109

#

870

azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths

1110

# azimuthAngle = arrayAOA[:,0]

871

A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))

1111

# zenithAngle = arrayAOA[:,1]

872

A = numpy.asmatrix(A)

1112

#

873

A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()

1113

# #Setting Error

874

velHor = numpy.dot(A1,velAux)

1114

# #Number 3: AOA not fesible

875

1115

# indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]

876

velEst[i,:] = numpy.squeeze(velHor)

1116

# error[indInvalid] = 3

877

return velEst

1117

# #Number 4: Large difference in AOAs obtained from different antenna baselines

1118

# indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]

1119

# error[indInvalid] = 4

1120

# return arrayAOA, error

1121

#

1122

# def __getDirectionCosines(self, arrayPhase, pairsList):

1123

#

1124

# #Initializing some variables

1125

# ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi

1126

# ang_aux = ang_aux.reshape(1,ang_aux.size)

1127

#

1128

# cosdir = numpy.zeros((arrayPhase.shape[0],2))

1129

# cosdir0 = numpy.zeros((arrayPhase.shape[0],2))

1130

#

1131

#

1132

# for i in range(2):

1133

# #First Estimation

1134

# phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]

1135

# #Dealias

1136

# indcsi = numpy.where(phi0_aux > numpy.pi)

1137

# phi0_aux[indcsi] -= 2*numpy.pi

1138

# indcsi = numpy.where(phi0_aux < -numpy.pi)

1139

# phi0_aux[indcsi] += 2*numpy.pi

1140

# #Direction Cosine 0

1141

# cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)

1142

#

1143

# #Most-Accurate Second Estimation

1144

# phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]

1145

# phi1_aux = phi1_aux.reshape(phi1_aux.size,1)

1146

# #Direction Cosine 1

1147

# cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)

1148

#

1149

# #Searching the correct Direction Cosine

1150

# cosdir0_aux = cosdir0[:,i]

1151

# cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)

1152

# #Minimum Distance

1153

# cosDiff = (cosdir1 - cosdir0_aux)**2

1154

# indcos = cosDiff.argmin(axis = 1)

1155

# #Saving Value obtained

1156

# cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]

1157

#

1158

# return cosdir0, cosdir

1159

#

1160

# def __calculateAOA(self, cosdir, azimuth):

1161

# cosdirX = cosdir[:,0]

1162

# cosdirY = cosdir[:,1]

1163

#

1164

# zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi

1165

# azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east

1166

# angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()

1167

#

1168

# return angles

1169

#

1170

# def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):

1171

#

1172

# Ramb = 375 #Ramb = c/(2*PRF)

1173

# Re = 6371 #Earth Radius

1174

# heights = numpy.zeros(Ranges.shape)

1175

#

1176

# R_aux = numpy.array([0,1,2])*Ramb

1177

# R_aux = R_aux.reshape(1,R_aux.size)

1178

#

1179

# Ranges = Ranges.reshape(Ranges.size,1)

1180

#

1181

# Ri = Ranges + R_aux

1182

# hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re

1183

#

1184

# #Check if there is a height between 70 and 110 km

1185

# h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)

1186

# ind_h = numpy.where(h_bool == 1)[0]

1187

#

1188

# hCorr = hi[ind_h, :]

1189

# ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))

1190

#

1191

# hCorr = hi[ind_hCorr]

1192

# heights[ind_h] = hCorr

1193

#

1194

# #Setting Error

1195

# #Number 13: Height unresolvable echo: not valid height within 70 to 110 km

1196

# #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km

1197

#

1198

# indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]

1199

# error[indInvalid2] = 14

1200

# indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]

1201

# error[indInvalid1] = 13

1202

#

1203

# return heights, error

1204

878

1205

def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):

879

def __getPhaseSlope(self, metArray, heightList, timeList):

880

meteorList = []

881

#utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2

882

#Putting back together the meteor matrix

883

utctime = metArray[:,0]

884

uniqueTime = numpy.unique(utctime)

1206

885

1207

'''

886

phaseDerThresh = 0.5

1208

Function GetMoments()

887

ippSeconds = timeList[1] - timeList[0]

888

sec = numpy.where(timeList>1)[0][0]

889

nPairs = metArray.shape[1] - 6

890

nHeights = len(heightList)

1209

891

1210

Input:

892

for t in uniqueTime:

1211

Output:

893

metArray1 = metArray[utctime==t,:]

1212

Variables modified:

894

# phaseDerThresh = numpy.pi/4 #reducir Phase thresh

1213

'''

895

tmet = metArray1[:,1].astype(int)

1214

if path != None:

896

hmet = metArray1[:,2].astype(int)

1215

sys.path.append(path)

897

1216

self.dataOut.library = importlib.import_module(file)

898

metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))

899

metPhase[:,:] = numpy.nan

900

metPhase[:,hmet,tmet] = metArray1[:,6:].T

901

902

#Delete short trails

903

metBool = ~numpy.isnan(metPhase[0,:,:])

904

heightVect = numpy.sum(metBool, axis = 1)

905

metBool[heightVect<sec,:] = False

906

metPhase[:,heightVect<sec,:] = numpy.nan

907

908

#Derivative

909

metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])

910

phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))

911

metPhase[phDerAux] = numpy.nan

912

913

#--------------------------METEOR DETECTION -----------------------------------------

914

indMet = numpy.where(numpy.any(metBool,axis=1))[0]

915

916

for p in numpy.arange(nPairs):

917

phase = metPhase[p,:,:]

918

phDer = metDer[p,:,:]

919

920

for h in indMet:

921

height = heightList[h]

922

phase1 = phase[h,:] #82

923

phDer1 = phDer[h,:]

924

925

phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap

926

927

indValid = numpy.where(~numpy.isnan(phase1))[0]

928

initMet = indValid[0]

929

endMet = 0

930

931

for i in range(len(indValid)-1):

932

933

#Time difference

934

inow = indValid[i]

935

inext = indValid[i+1]

936

idiff = inext - inow

937

#Phase difference

938

phDiff = numpy.abs(phase1[inext] - phase1[inow])

939

940

if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor

941

sizeTrail = inow - initMet + 1

942

if sizeTrail>3*sec: #Too short meteors

943

x = numpy.arange(initMet,inow+1)*ippSeconds

944

y = phase1[initMet:inow+1]

945

ynnan = ~numpy.isnan(y)

946

x = x[ynnan]

947

y = y[ynnan]

948

slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)

949

ylin = x*slope + intercept

950

rsq = r_value**2

951

if rsq > 0.5:

952

vel = slope#*height*1000/(k*d)

953

estAux = numpy.array([utctime,p,height, vel, rsq])

954

meteorList.append(estAux)

955

initMet = inext

956

metArray2 = numpy.array(meteorList)

1217

957

1218

#To be inserted as a parameter

958

return metArray2

1219

groupArray = numpy.array(groupList)

959

1220

# groupArray = numpy.array([[0,1],[2,3]])

960

def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):

1221

self.dataOut.groupList = groupArray

1222

961

1223

nGroups = groupArray.shape[0]

962

azimuth1 = numpy.zeros(len(pairslist))

1224

nChannels = self.dataIn.nChannels

963

dist = numpy.zeros(len(pairslist))

1225

nHeights=self.dataIn.heightList.size

1226

964

1227

#Parameters Array

965

for i in range(len(rx_location)):

1228

self.dataOut.data_param = None

966

ch0 = pairslist[i][0]

967

ch1 = pairslist[i][1]

968

969

diffX = rx_location[ch0][0] - rx_location[ch1][0]

970

diffY = rx_location[ch0][1] - rx_location[ch1][1]

971

azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi

972

dist[i] = numpy.sqrt(diffX**2 + diffY**2)

1229

973

1230

#Set constants

974

azimuth1 -= azimuth0

1231

constants = self.dataOut.library.setConstants(self.dataIn)

975

return azimuth1, dist

1232

self.dataOut.constants = constants

976

1233

M = self.dataIn.normFactor

977

def techniqueNSM_DBS(self, **kwargs):

1234

N = self.dataIn.nFFTPoints

978

metArray = kwargs['metArray']

1235

ippSeconds = self.dataIn.ippSeconds

979

heightList = kwargs['heightList']

1236

K = self.dataIn.nIncohInt

980

timeList = kwargs['timeList']

1237

pairsArray = numpy.array(self.dataIn.pairsList)

981

zenithList = kwargs['zenithList']

982

nChan = numpy.max(cmet) + 1

983

nHeights = len(heightList)

1238

984

1239

#List of possible combinations

985

utctime = metArray[:,0]

1240

listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)

986

cmet = metArray[:,1]

1241

indCross = numpy.zeros(len(list(listComb)), dtype = 'int')

987

hmet = metArray1[:,3].astype(int)

988

h1met = heightList[hmet]*zenithList[cmet]

989

vmet = metArray1[:,5]

1242

990

1243

if getSNR:

991

for i in range(nHeights - 1):

1244

listChannels = groupArray.reshape((groupArray.size))

992

hmin = heightList[i]

1245

listChannels.sort()

993

hmax = heightList[i + 1]

1246

noise = self.dataIn.getNoise()

994

1247

self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])

995

vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]

996

997

998

999

return data_output

1000

1001

def run(self, dataOut, technique, **kwargs):

1248

1002

1249

for i in range(nGroups):

1003

param = dataOut.data_param

1250

coord = groupArray[i,:]

1004

if dataOut.abscissaList != None:

1005

absc = dataOut.abscissaList[:-1]

1006

noise = dataOut.noise

1007

heightList = dataOut.heightList

1008

SNR = dataOut.data_SNR

1009

1010

if technique == 'DBS':

1011

# if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):

1012

# theta_x = numpy.array(kwargs['dirCosx'])

1013

# theta_y = numpy.array(kwargs['dirCosy'])

1014

# else:

1015

# elev = numpy.array(kwargs['elevation'])

1016

# azim = numpy.array(kwargs['azimuth'])

1017

# theta_x, theta_y = self.__calculateCosDir(elev, azim)

1018

# azimuth = kwargs['correctAzimuth']

1019

# if kwargs.has_key('horizontalOnly'):

1020

# horizontalOnly = kwargs['horizontalOnly']

1021

# else: horizontalOnly = False

1022

# if kwargs.has_key('correctFactor'):

1023

# correctFactor = kwargs['correctFactor']

1024

# else: correctFactor = 1

1025

# if kwargs.has_key('channelList'):

1026

# channelList = kwargs['channelList']

1027

# if len(channelList) == 2:

1028

# horizontalOnly = True

1029

# arrayChannel = numpy.array(channelList)

1030

# param = param[arrayChannel,:,:]

1031

# theta_x = theta_x[arrayChannel]

1032

# theta_y = theta_y[arrayChannel]

1033

1034

velRadial0 = param[:,1,:] #Radial velocity

1035

# dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function

1036

dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, heightList, SNR, kwargs) #DBS Function

1037

dataOut.utctimeInit = dataOut.utctime

1038

dataOut.outputInterval = dataOut.paramInterval

1251

1039

1252

#Input data array

1040

elif technique == 'SA':

1253

data = self.dataIn.data_spc[coord,:,:]/(M*N)

1041

1254

data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))

1042

#Parameters

1043

position_x = kwargs['positionX']

1044

position_y = kwargs['positionY']

1045

azimuth = kwargs['azimuth']

1255

1046

1256

#Cross Spectra data array for Covariance Matrixes

1047

if kwargs.has_key('crosspairsList'):

1257

ind = 0

1048

pairs = kwargs['crosspairsList']

1258

for pairs in listComb:

1049

else:

1259

pairsSel = numpy.array([coord[x],coord[y]])

1050

pairs = None

1260

indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])

1051

1261

ind += 1

1052

if kwargs.has_key('correctFactor'):

1262

dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)

1053

correctFactor = kwargs['correctFactor']

1263

dataCross = dataCross**2/K

1054

else:

1055

correctFactor = 1

1056

1057

tau = dataOut.data_param

1058

_lambda = dataOut.C/dataOut.frequency

1059

pairsList = dataOut.groupList

1060

nChannels = dataOut.nChannels

1061

1062

dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)

1063

dataOut.utctimeInit = dataOut.utctime

1064

dataOut.outputInterval = dataOut.timeInterval

1264

1065

1265

for h in range(nHeights):

1066

elif technique == 'Meteors':

1266

# print self.dataOut.heightList[h]

1067

dataOut.flagNoData = True

1068

self.__dataReady = False

1069

1070

if kwargs.has_key('nHours'):

1071

nHours = kwargs['nHours']

1072

else:

1073

nHours = 1

1267

1074

1268

#Input

1075

if kwargs.has_key('meteorsPerBin'):

1269

d = data[:,h]

1076

meteorThresh = kwargs['meteorsPerBin']

1077

else:

1078

meteorThresh = 6

1079

1080

if kwargs.has_key('hmin'):

1081

hmin = kwargs['hmin']

1082

else: hmin = 70

1083

if kwargs.has_key('hmax'):

1084

hmax = kwargs['hmax']

1085

else: hmax = 110

1086

1087

dataOut.outputInterval = nHours*3600

1088

1089

if self.__isConfig == False:

1090

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

1091

#Get Initial LTC time

1092

self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

1093

self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

1270

1094

1271

#Covariance Matrix

1095

self.__isConfig = True

1272

D = numpy.diag(d**2/K)

1096

1273

ind = 0

1097

if self.__buffer == None:

1274

for pairs in listComb:

1098

self.__buffer = dataOut.data_param

1275

#Coordinates in Covariance Matrix

1099

self.__firstdata = copy.copy(dataOut)

1276

x = pairs[0]

1277

y = pairs[1]

1278

#Channel Index

1279

S12 = dataCross[ind,:,h]

1280

D12 = numpy.diag(S12)

1281

#Completing Covariance Matrix with Cross Spectras

1282

D[x*N:(x+1)*N,y*N:(y+1)*N] = D12

1283

D[y*N:(y+1)*N,x*N:(x+1)*N] = D12

1284

ind += 1

1285

Dinv=numpy.linalg.inv(D)

1286

L=numpy.linalg.cholesky(Dinv)

1287

LT=L.T

1288

1100

1289

dp = numpy.dot(LT,d)

1101

else:

1102

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

1103

1104

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

1105

1106

if self.__dataReady:

1107

dataOut.utctimeInit = self.__initime

1290

1108

1291

#Initial values

1109

self.__initime += dataOut.outputInterval #to erase time offset

1292

data_spc = self.dataIn.data_spc[coord,:,h]

1293

1110

1294

if (h>0)and(error1[3]<5):

1111

dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)

1295

p0 = self.dataOut.data_param[i,:,h-1]

1112

dataOut.flagNoData = False

1296

~~else~~:

1113

self.__buffer = None

1297

p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))

1298

1114

1299

try:

1115

elif technique == 'Meteors1':

1300

#Least Squares

1116

dataOut.flagNoData = True

1301

minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)

1117

self.__dataReady = False

1302

# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))

1118

1303

#Chi square error

1119

if kwargs.has_key('nMins'):

1304

error0 = numpy.sum(infodict['fvec']**2)/(2*N)

1120

nMins = kwargs['nMins']

1305

#Error with Jacobian

1121

else: nMins = 20

1306

error1 = self.dataOut.library.errorFunction(minp,constants,LT)

1122

if kwargs.has_key('rx_location'):

1307

except:

1123

rx_location = kwargs['rx_location']

1308

minp = p0*numpy.nan

1124

else: rx_location = [(0,1),(1,1),(1,0)]

1309

error0 = numpy.nan

1125

if kwargs.has_key('azimuth'):

1310

error1 = p0*numpy.nan

1126

azimuth = kwargs['azimuth']

1311

1127

else: azimuth = 51

1312

#Save

1128

if kwargs.has_key('dfactor'):

1313

if self.dataOut.data_param == None:

1129

dfactor = kwargs['dfactor']

1314

self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan

1130

if kwargs.has_key('mode'):

1315

self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan

1131

mode = kwargs['mode']

1132

else: mode = 'SA'

1133

1134

#Borrar luego esto

1135

if dataOut.groupList == None:

1136

dataOut.groupList = [(0,1),(0,2),(1,2)]

1137

groupList = dataOut.groupList

1138

C = 3e8

1139

freq = 50e6

1140

lamb = C/freq

1141

k = 2*numpy.pi/lamb

1142

1143

timeList = dataOut.abscissaList

1144

heightList = dataOut.heightList

1145

1146

if self.__isConfig == False:

1147

dataOut.outputInterval = nMins*60

1148

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

1149

#Get Initial LTC time

1150

initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

1151

minuteAux = initime.minute

1152

minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)

1153

self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

1154

1155

self.__isConfig = True

1316

1156

1317

self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))

1157

if self.__buffer == None:

1318

self.dataOut.data_param[i~~,:,~~h] = ~~minp~~

1158

self.__buffer = dataOut.data_param

1159

self.__firstdata = copy.copy(dataOut)

1160

1161

else:

1162

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

1163

1164

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

1165

1166

if self.__dataReady:

1167

dataOut.utctimeInit = self.__initime

1168

self.__initime += dataOut.outputInterval #to erase time offset

1169

1170

metArray = self.__buffer

1171

if mode == 'SA':

1172

dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)

1173

elif mode == 'DBS':

1174

dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)

1175

dataOut.data_output = dataOut.data_output.T

1176

dataOut.flagNoData = False

1177

self.__buffer = None

1178

1179

return

1180

1181

class EWDriftsEstimation(Operation):

1182

1183

def __init__(self):

1184

Operation.__init__(self)

1185

1186

def __correctValues(self, heiRang, phi, velRadial, SNR):

1187

listPhi = phi.tolist()

1188

maxid = listPhi.index(max(listPhi))

1189

minid = listPhi.index(min(listPhi))

1190

1191

rango = range(len(phi))

1192

# rango = numpy.delete(rango,maxid)

1193

1194

heiRang1 = heiRang*math.cos(phi[maxid])

1195

heiRangAux = heiRang*math.cos(phi[minid])

1196

indOut = (heiRang1 < heiRangAux[0]).nonzero()

1197

heiRang1 = numpy.delete(heiRang1,indOut)

1198

1199

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

1200

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

1201

1202

for i in rango:

1203

x = heiRang*math.cos(phi[i])

1204

y1 = velRadial[i,:]

1205

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

1206

1207

x1 = heiRang1

1208

y11 = f1(x1)

1209

1210

y2 = SNR[i,:]

1211

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

1212

y21 = f2(x1)

1213

1214

velRadial1[i,:] = y11

1215

SNR1[i,:] = y21

1216

1217

return heiRang1, velRadial1, SNR1

1218

1219

def run(self, dataOut, zenith, zenithCorrection):

1220

heiRang = dataOut.heightList

1221

velRadial = dataOut.data_param[:,3,:]

1222

SNR = dataOut.data_SNR

1223

1224

zenith = numpy.array(zenith)

1225

zenith -= zenithCorrection

1226

zenith *= numpy.pi/180

1227

1228

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)

1229

1230

alp = zenith[0]

1231

bet = zenith[1]

1232

1233

w_w = velRadial1[0,:]

1234

w_e = velRadial1[1,:]

1235

1236

w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))

1237

u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))

1238

1239

winds = numpy.vstack((u,w))

1240

1241

dataOut.heightList = heiRang1

1242

dataOut.data_output = winds

1243

dataOut.data_SNR = SNR1

1244

1245

dataOut.utctimeInit = dataOut.utctime

1246

dataOut.outputInterval = dataOut.timeInterval

1319

return

1247

return

1320

1321

def __residFunction(self, p, dp, LT, constants):

1322

1248

1323

fm = self.dataOut.library.modelFunction(p, constants)

1249

#--------------- Non Specular Meteor ----------------

1324

fmp=numpy.dot(LT,fm)

1325

1326

return dp-fmp

1327

1250

1328

def __getSNR(self, z, noise):

1251

class NonSpecularMeteorDetection(Operation):

1329

1252

1330

avg = numpy.average(z, axis=1)

1253

def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):

1331

SNR = (avg.T-noise)/noise

1332

SNR = SNR.T

1333

return SNR

1334

1335

def __chisq(p,chindex,hindex):

1336

#similar to Resid but calculates CHI**2

1337

[LT,d,fm]=setupLTdfm(p,chindex,hindex)

1338

dp=numpy.dot(LT,d)

1339

fmp=numpy.dot(LT,fm)

1340

chisq=numpy.dot((dp-fmp).T,(dp-fmp))

1341

return chisq

1342

1343

def NonSpecularMeteorDetection(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):

1344

data_acf = self.dataOut.data_pre[0]

1254

data_acf = self.dataOut.data_pre[0]

1345

data_ccf = self.dataOut.data_pre[1]

1255

data_ccf = self.dataOut.data_pre[1]

1346

1256

1485

coordMet = numpy.where(boolMetFin)

1395

coordMet = numpy.where(boolMetFin)

1486

1396

1487

cmet = coordMet[0]

1397

cmet = coordMet[0]

1488

tmet = coordMet[1]

1398

tmet = coordMet[1]

1489

hmet = coordMet[2]

1399

hmet = coordMet[2]

1490

1491

data_param = numpy.zeros((tmet.size, 7))

1492

data_param[:,0] = utctime

1493

data_param[:,1] = cmet

1494

data_param[:,2] = tmet

1495

data_param[:,3] = hmet

1496

data_param[:,4] = SNR[cmet,tmet,hmet].T

1497

data_param[:,5] = velRad[cmet,tmet,hmet].T

1498

data_param[:,6] = spcWidth[cmet,tmet,hmet].T

1499

1500

# self.dataOut.data_param = data_int

1501

if len(data_param) == 0:

1502

self.dataOut.flagNoData = True

1503

else:

1504

self.dataOut.data_param = data_param

1505

1506

def __erase_small(self, binArray, threshX, threshY):

1507

labarray, numfeat = ndimage.measurements.label(binArray)

1508

binArray1 = numpy.copy(binArray)

1509

1510

for i in range(1,numfeat + 1):

1511

auxBin = (labarray==i)

1512

auxSize = auxBin.sum()

1513

1514

x,y = numpy.where(auxBin)

1515

widthX = x.max() - x.min()

1516

widthY = y.max() - y.min()

1517

1518

#width X: 3 seg -> 12.5*3

1519

#width Y:

1520

1521

if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):

1522

binArray1[auxBin] = False

1523

1524

return binArray1

1525

1526

def WeirdEcho(self):

1527

# data_pre = self.dataOut.data_pre

1528

# nHeights = self.dataOut.nHeights

1529

# # nProfiles = self.dataOut.data_pre.shape[1]

1530

# # data_param = numpy.zeros((len(pairslist),nProfiles,nHeights))

1531

# data_param = numpy.zeros((len(pairslist),nHeights))

1532

#

1533

# for i in range(len(pairslist)):

1534

# chan0 = data_pre[pairslist[i][0],:]

1535

# chan1 = data_pre[pairslist[i][1],:]

1536

# #calcular correlacion cruzada

1537

# #magnitud es coherencia

1538

# #fase es dif fase

1539

# correl = chan0*numpy.conj(chan1)

1540

# coherence = numpy.abs(correl)/(numpy.abs(chan0)*numpy.abs(chan1))

1541

# phase = numpy.angle(correl)

1542

# # data_param[2*i,:,:] = coherence

1543

# data_param[i,:] = phase

1544

#

1545

# self.dataOut.data_param = data_param

1546

# self.dataOut.groupList = pairslist

1547

data_cspc = self.dataOut.data_pre[1]

1548

ccf = numpy.average(data_cspc,axis=1)

1549

phases = numpy.angle(ccf).T

1550

1551

meteorOps = MeteorOperations()

1552

pairsList = ((0,1),(2,3))

1553

jph = numpy.array([0,0,0,0])

1554

AOAthresh = numpy.pi/8

1555

azimuth = 45

1556

error = numpy.zeros((phases.shape[0],1))

1557

AOA,error = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)

1558

self.dataOut.data_param = AOA.T

1559

1560

class WindProfiler(Operation):

1561

1562

__isConfig = False

1563

1564

__initime = None

1565

__lastdatatime = None

1566

__integrationtime = None

1567

1568

__buffer = None

1569

1570

__dataReady = False

1571

1572

__firstdata = None

1573

1574

n = None

1575

1576

def __init__(self):

1577

Operation.__init__(self)

1578

1579

def __calculateCosDir(self, elev, azim):

1580

zen = (90 - elev)*numpy.pi/180

1581

azim = azim*numpy.pi/180

1582

cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))

1583

cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)

1584

1585

signX = numpy.sign(numpy.cos(azim))

1586

signY = numpy.sign(numpy.sin(azim))

1587

1588

cosDirX = numpy.copysign(cosDirX, signX)

1589

cosDirY = numpy.copysign(cosDirY, signY)

1590

return cosDirX, cosDirY

1591

1592

def __calculateAngles(self, theta_x, theta_y, azimuth):

1593

1594

dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)

1595

zenith_arr = numpy.arccos(dir_cosw)

1596

azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180

1597

1598

dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)

1599

dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)

1600

1601

return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw

1602

1603

def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):

1604

1605

#

1606

if horOnly:

1607

A = numpy.c_[dir_cosu,dir_cosv]

1608

else:

1609

A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]

1610

A = numpy.asmatrix(A)

1611

A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()

1612

1613

return A1

1614

1615

def __correctValues(self, heiRang, phi, velRadial, SNR):

1616

listPhi = phi.tolist()

1617

maxid = listPhi.index(max(listPhi))

1618

minid = listPhi.index(min(listPhi))

1619

1620

rango = range(len(phi))

1621

# rango = numpy.delete(rango,maxid)

1622

1623

heiRang1 = heiRang*math.cos(phi[maxid])

1624

heiRangAux = heiRang*math.cos(phi[minid])

1625

indOut = (heiRang1 < heiRangAux[0]).nonzero()

1626

heiRang1 = numpy.delete(heiRang1,indOut)

1627

1628

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

1629

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

1630

1631

for i in rango:

1632

x = heiRang*math.cos(phi[i])

1633

y1 = velRadial[i,:]

1634

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

1635

1636

x1 = heiRang1

1637

y11 = f1(x1)

1638

1639

y2 = SNR[i,:]

1640

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

1641

y21 = f2(x1)

1642

1643

velRadial1[i,:] = y11

1644

SNR1[i,:] = y21

1645

1646

return heiRang1, velRadial1, SNR1

1647

1648

def __calculateVelUVW(self, A, velRadial):

1649

1650

#Operacion Matricial

1651

# velUVW = numpy.zeros((velRadial.shape[1],3))

1652

# for ind in range(velRadial.shape[1]):

1653

# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])

1654

# velUVW = velUVW.transpose()

1655

velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))

1656

velUVW[:,:] = numpy.dot(A,velRadial)

1657

1658

1659

return velUVW

1660

1661

def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):

1662

"""

1663

Function that implements Doppler Beam Swinging (DBS) technique.

1664

1665

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

1666

Direction correction (if necessary), Ranges and SNR

1667

1668

Output: Winds estimation (Zonal, Meridional and Vertical)

1669

1670

Parameters affected: Winds, height range, SNR

1671

"""

1672

azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)

1673

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)

1674

A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)

1675

1676

#Calculo de Componentes de la velocidad con DBS

1677

winds = self.__calculateVelUVW(A,velRadial1)

1678

1679

return winds, heiRang1, SNR1

1680

1681

def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):

1682

1683

posx = numpy.asarray(posx)

1684

posy = numpy.asarray(posy)

1685

1686

#Rotacion Inversa para alinear con el azimuth

1687

if azimuth!= None:

1688

azimuth = azimuth*math.pi/180

1689

posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)

1690

posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)

1691

else:

1692

posx1 = posx

1693

posy1 = posy

1694

1695

#Calculo de Distancias

1696

distx = numpy.zeros(pairsCrossCorr.size)

1697

disty = numpy.zeros(pairsCrossCorr.size)

1698

dist = numpy.zeros(pairsCrossCorr.size)

1699

ang = numpy.zeros(pairsCrossCorr.size)

1700

1701

for i in range(pairsCrossCorr.size):

1702

distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]

1703

disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]

1704

dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)

1705

ang[i] = numpy.arctan2(disty[i],distx[i])

1706

#Calculo de Matrices

1707

nPairs = len(pairs)

1708

ang1 = numpy.zeros((nPairs, 2, 1))

1709

dist1 = numpy.zeros((nPairs, 2, 1))

1710

1711

for j in range(nPairs):

1712

dist1[j,0,0] = dist[pairs[j][0]]

1713

dist1[j,1,0] = dist[pairs[j][1]]

1714

ang1[j,0,0] = ang[pairs[j][0]]

1715

ang1[j,1,0] = ang[pairs[j][1]]

1716

1400

1717

return distx,disty, dist1,ang1

1401

data_param = numpy.zeros((tmet.size, 7))

1718

1402

data_param[:,0] = utctime

1719

def __calculateVelVer(self, phase, lagTRange, _lambda):

1403

data_param[:,1] = cmet

1404

data_param[:,2] = tmet

1405

data_param[:,3] = hmet

1406

data_param[:,4] = SNR[cmet,tmet,hmet].T

1407

data_param[:,5] = velRad[cmet,tmet,hmet].T

1408

data_param[:,6] = spcWidth[cmet,tmet,hmet].T

1409

1410

# self.dataOut.data_param = data_int

1411

if len(data_param) == 0:

1412

self.dataOut.flagNoData = True

1413

else:

1414

self.dataOut.data_param = data_param

1720

1415

1721

Ts = lagTRange[1] - lagTRange[0]

1416

def __erase_small(self, binArray, threshX, threshY):

1722

velW = -_lambda*phase/(4*math.pi*Ts)

1417

labarray, numfeat = ndimage.measurements.label(binArray)

1723

1418

binArray1 = numpy.copy(binArray)

1724

return velW

1725

1726

def __calculateVelHorDir(self, dist, tau1, tau2, ang):

1727

nPairs = tau1.shape[0]

1728

vel = numpy.zeros((nPairs,3,tau1.shape[2]))

1729

1419

1730

angCos = numpy.cos(ang)

1420

for i in range(1,numfeat + 1):

1731

angSin = numpy.sin(ang)

1421

auxBin = (labarray==i)

1422

auxSize = auxBin.sum()

1423

1424

x,y = numpy.where(auxBin)

1425

widthX = x.max() - x.min()

1426

widthY = y.max() - y.min()

1427

1428

#width X: 3 seg -> 12.5*3

1429

#width Y:

1430

1431

if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):

1432

binArray1[auxBin] = False

1433

1434

return binArray1

1435

1436

#--------------- Specular Meteor ----------------

1437

1438

class MeteorDetection(Operation):

1439

'''

1440

Function DetectMeteors()

1441

Project developed with paper:

1442

HOLDSWORTH ET AL. 2004

1443

1444

Input:

1445

self.dataOut.data_pre

1732

1446

1733

vel0 = dist*tau1/(2*tau2**2)

1447

centerReceiverIndex: From the channels, which is the center receiver

1734

vel[:,0,:] = (vel0*angCos).sum(axis = 1)

1448

1735

vel[:,1,:] = (vel0*angSin).sum(axis = 1)

1449

hei_ref: Height reference for the Beacon signal extraction

1450

tauindex:

1451

predefinedPhaseShifts: Predefined phase offset for the voltge signals

1452

1453

cohDetection: Whether to user Coherent detection or not

1454

cohDet_timeStep: Coherent Detection calculation time step

1455

cohDet_thresh: Coherent Detection phase threshold to correct phases

1456

1457

noise_timeStep: Noise calculation time step

1458

noise_multiple: Noise multiple to define signal threshold

1459

1460

multDet_timeLimit: Multiple Detection Removal time limit in seconds

1461

multDet_rangeLimit: Multiple Detection Removal range limit in km

1462

1463

phaseThresh: Maximum phase difference between receiver to be consider a meteor

1464

SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor

1465

1466

hmin: Minimum Height of the meteor to use it in the further wind estimations

1467

hmax: Maximum Height of the meteor to use it in the further wind estimations

1468

azimuth: Azimuth angle correction

1469

1470

Affected:

1471

self.dataOut.data_param

1736

1472

1737

ind = numpy.where(numpy.isinf(vel))

1473

Rejection Criteria (Errors):

1738

vel[ind] = numpy.nan

1474

0: No error; analysis OK

1739

1475

1: SNR < SNR threshold

1740

return vel

1476

2: angle of arrival (AOA) ambiguously determined

1741

1477

3: AOA estimate not feasible

1742

def __getPairsAutoCorr(self, pairsList, nChannels):

1478

4: Large difference in AOAs obtained from different antenna baselines

1743

1479

5: echo at start or end of time series

1744

pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

1480

6: echo less than 5 examples long; too short for analysis

1481

7: echo rise exceeds 0.3s

1482

8: echo decay time less than twice rise time

1483

9: large power level before echo

1484

10: large power level after echo

1485

11: poor fit to amplitude for estimation of decay time

1486

12: poor fit to CCF phase variation for estimation of radial drift velocity

1487

13: height unresolvable echo: not valid height within 70 to 110 km

1488

14: height ambiguous echo: more then one possible height within 70 to 110 km

1489

15: radial drift velocity or projected horizontal velocity exceeds 200 m/s

1490

16: oscilatory echo, indicating event most likely not an underdense echo

1745

1491

1746

for l in range(len(pairsList)):

1492

17: phase difference in meteor Reestimation

1747

firstChannel = pairsList[l][0]

1748

secondChannel = pairsList[l][1]

1749

1750

#Obteniendo pares de Autocorrelacion

1751

if firstChannel == secondChannel:

1752

pairsAutoCorr[firstChannel] = int(l)

1753

1754

pairsAutoCorr = pairsAutoCorr.astype(int)

1755

1493

1756

pairsCrossCorr = range(len(pairsList))

1494

Data Storage:

1757

pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

1495

Meteors for Wind Estimation (8):

1496

Utc Time | Range Height

1497

Azimuth Zenith errorCosDir

1498

VelRad errorVelRad

1499

Phase0 Phase1 Phase2 Phase3

1500

TypeError

1758

1501

1759

return pairsAutoCorr, pairsCrossCorr

1502

'''

1760

1503

1761

def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):

1504

def run(self, dataOut, hei_ref = None, tauindex = 0,

1762

"""

1505

phaseOffsets = None,

1763

Function that implements Spaced Antenna (SA) technique.

1506

cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,

1764

1507

noise_timeStep = 4, noise_multiple = 4,

1765

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

1508

multDet_timeLimit = 1, multDet_rangeLimit = 3,

1766

Direction correction (if necessary), Ranges and SNR

1509

phaseThresh = 20, SNRThresh = 5,

1510

hmin = 50, hmax=150, azimuth = 0,

1511

channelPositions = None) :

1767

1512

1768

Output: Winds estimation (Zonal, Meridional and Vertical)

1513

1514

#Getting Pairslist

1515

if channelPositions == None:

1516

# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T

1517

channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella

1518

meteorOps = MeteorOperations()

1519

pairslist0, distances = meteorOps.getPhasePairs(channelPositions)

1769

1520

1770

Parameters affected: Winds

1521

#Get Beacon signal

1771

"""

1522

newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])

1772

#Cross Correlation pairs obtained

1773

pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

1774

pairsArray = numpy.array(pairsList)[pairsCrossCorr]

1775

pairsSelArray = numpy.array(pairsSelected)

1776

pairs = []

1777

1523

1778

#Wind estimation pairs obtained

1524

if hei_ref != None:

1779

for i in range(pairsSelArray.shape[0]/2):

1525

newheis = numpy.where(self.dataOut.heightList>hei_ref)

1780

ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]

1781

ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]

1782

pairs.append((ind1,ind2))

1783

1526

1784

indtau = tau.shape[0]/2

1527

heiRang = dataOut.getHeiRange()

1785

tau1 = tau[:indtau,:]

1528

1786

tau2 = tau[indtau:-1,:]

1529

#****************REMOVING HARDWARE PHASE DIFFERENCES***************

1787

tau1 = tau1[pairs,:]

1530

# see if the user put in pre defined phase shifts

1788

tau2 = tau2[pairs,:]

1531

voltsPShift = self.dataOut.data_pre.copy()

1789

phase1 = tau[-1,:]

1790

1532

1791

#---------------------------------------------------------------------

1533

# if predefinedPhaseShifts != None:

1792

#Metodo Directo

1534

# hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180

1793

distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)

1535

#

1794

winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)

1536

# # elif beaconPhaseShifts:

1795

winds = stats.nanmean(winds, axis=0)

1537

# # #get hardware phase shifts using beacon signal

1796

#---------------------------------------------------------------------

1538

# # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)

1797

#Metodo General

1539

# # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)

1798

# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)

1540

#

1799

# #Calculo Coeficientes de Funcion de Correlacion

1541

# else:

1800

# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)

1542

# hardwarePhaseShifts = numpy.zeros(5)

1801

# #Calculo de Velocidades

1543

#

1802

# winds = self.calculateVelUV(F,G,A,B,H)

1544

# voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')

1545

# for i in range(self.dataOut.data_pre.shape[0]):

1546

# voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])

1803

1547

1804

#---------------------------------------------------------------------

1548

#******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********

1805

winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)

1806

winds = correctFactor*winds

1807

return winds

1808

1809

def __checkTime(self, currentTime, paramInterval, outputInterval):

1810

1811

dataTime = currentTime + paramInterval

1812

deltaTime = dataTime - self.__initime

1813

1814

if deltaTime >= outputInterval or deltaTime < 0:

1815

self.__dataReady = True

1816

return

1817

1549

1818

def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):

1550

#Remove DC

1819

'''

1551

voltsDC = numpy.mean(voltsPShift,1)

1820

Function that implements winds estimation technique with detected meteors.

1552

voltsDC = numpy.mean(voltsDC,1)

1821

1553

for i in range(voltsDC.shape[0]):

1822

Input: Detected meteors, Minimum meteor quantity to wind estimation

1554

voltsPShift[i] = voltsPShift[i] - voltsDC[i]

1823

1555

1824

Output: Winds estimation (Zonal and Meridional)

1556

#Don't considerate last heights, theyre used to calculate Hardware Phase Shift

1825

1557

voltsPShift = voltsPShift[:,:,:newheis[0][0]]

1826

Parameters affected: Winds

1827

'''

1828

# print arrayMeteor.shape

1829

#Settings

1830

nInt = (heightMax - heightMin)/2

1831

# print nInt

1832

nInt = int(nInt)

1833

# print nInt

1834

winds = numpy.zeros((2,nInt))*numpy.nan

1835

1836

#Filter errors

1837

error = numpy.where(arrayMeteor[:,-1] == 0)[0]

1838

finalMeteor = arrayMeteor[error,:]

1839

1558

1840

#Meteor Histogram

1559

#************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********

1841

finalHeights = finalMeteor[:,2]

1560

#Coherent Detection

1842

hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))

1561

if cohDetection:

1843

nMeteorsPerI = hist[0]

1562

#use coherent detection to get the net power

1844

heightPerI = hist[1]

1563

cohDet_thresh = cohDet_thresh*numpy.pi/180

1564

voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh)

1845

1565

1846

#~~Sort of meteors~~

1566

#Non-coherent detection!

1847

indSort = finalHeights.argsort()

1567

powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)

1848

finalMeteor2 = finalMeteor[indSort,:]

1568

#********** END OF COH/NON-COH POWER CALCULATION**********************

1569

1570

#********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************

1571

#Get noise

1572

noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)

1573

# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)

1574

#Get signal threshold

1575

signalThresh = noise_multiple*noise

1576

#Meteor echoes detection

1577

listMeteors = self.__findMeteors(powerNet, signalThresh)

1578

#******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********

1849

1579

1850

# Calculating winds

1580

#************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************

1851

ind1 = 0

1581

#Parameters

1852

ind2 = 0

1582

heiRange = self.dataOut.getHeiRange()

1583

rangeInterval = heiRange[1] - heiRange[0]

1584

rangeLimit = multDet_rangeLimit/rangeInterval

1585

timeLimit = multDet_timeLimit/self.dataOut.timeInterval

1586

#Multiple detection removals

1587

listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)

1588

#************ END OF REMOVE MULTIPLE DETECTIONS **********************

1853

1589

1854

for i in range(nInt):

1590

#********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************

1855

nMet = nMeteorsPerI[i]

1591

#Parameters

1856

ind1 = ind2

1592

phaseThresh = phaseThresh*numpy.pi/180

1857

ind2 = ind1 + nMet

1593

thresh = [phaseThresh, noise_multiple, SNRThresh]

1594

#Meteor reestimation (Errors N 1, 6, 12, 17)

1595

listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)

1596

# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)

1597

#Estimation of decay times (Errors N 7, 8, 11)

1598

listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)

1599

#******************* END OF METEOR REESTIMATION *******************

1600

1601

#********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************

1602

#Calculating Radial Velocity (Error N 15)

1603

radialStdThresh = 10

1604

listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval)

1605

1606

if len(listMeteors4) > 0:

1607

#Setting New Array

1608

date = self.dataOut.utctime

1609

arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)

1858

1610

1859

meteorAux = finalMeteor2[ind1:ind2,:]

1611

#Correcting phase offset

1612

if phaseOffsets != None:

1613

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

1614

arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

1860

1615

1861

if meteorAux.shape[0] >= meteorThresh:

1616

#Second Pairslist

1862

vel = meteorAux[:, 6]

1617

pairsList = []

1863

zen = meteorAux[:, 4]*numpy.pi/180

1618

pairx = (0,1)

1864

azim = meteorAux[:, 3]*numpy.pi/180

1619

pairy = (2,3)

1865

1620

pairsList.append(pairx)

1866

n = numpy.cos(zen)

1621

pairsList.append(pairy)

1867

# m = (1 - n**2)/(1 - numpy.tan(azim)**2)

1622

1868

# l = m*numpy.tan(azim)

1623

jph = numpy.array([0,0,0,0])

1869

l = numpy.sin(zen)*numpy.sin(azim)

1624

h = (hmin,hmax)

1870

m = numpy.sin(zen)*numpy.cos(azim)

1625

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

1871

1626

1872

A = numpy.vstack((l, m)).transpose()

1627

# #Calculate AOA (Error N 3, 4)

1873

A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())

1628

# #JONES ET AL. 1998

1874

windsAux = numpy.dot(A1, vel)

1629

# error = arrayParameters[:,-1]

1875

1630

# AOAthresh = numpy.pi/8

1876

winds[0,i] = windsAux[0]

1631

# phases = -arrayParameters[:,9:13]

1877

winds[1,i] = windsAux[1]

1632

# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)

1633

#

1634

# #Calculate Heights (Error N 13 and 14)

1635

# error = arrayParameters[:,-1]

1636

# Ranges = arrayParameters[:,2]

1637

# zenith = arrayParameters[:,5]

1638

# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)

1639

# error = arrayParameters[:,-1]

1640

#********************* END OF PARAMETERS CALCULATION **************************

1878

1641

1879

return winds, heightPerI[:-1]

1642

#***************************+ PASS DATA TO NEXT STEP **********************

1880

1643

# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))

1881

def techniqueNSM_SA(self, **kwargs):

1644

self.dataOut.data_param = arrayParameters

1882

metArray = kwargs['metArray']

1645

1883

heightList = kwargs['heightList']

1646

if arrayParameters == None:

1884

timeList = kwargs['timeList']

1647

self.dataOut.flagNoData = True

1648

else:

1649

self.dataOut.flagNoData = True

1650

1651

return

1885

1652

1886

rx_location = kwargs['rx_location']

1653

def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):

1887

groupList = kwargs['groupList']

1654

1888

azimuth = kwargs['azimuth']

1655

minIndex = min(newheis[0])

1889

dfactor = kwargs['dfactor']

1656

maxIndex = max(newheis[0])

1890

k = kwargs['k']

1891

1657

1892

azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)

1658

voltage = voltage0[:,:,minIndex:maxIndex+1]

1893

d = dist*dfactor

1659

nLength = voltage.shape[1]/n

1894

#Phase calculation

1660

nMin = 0

1895

metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)

1661

nMax = 0

1662

phaseOffset = numpy.zeros((len(pairslist),n))

1896

1663

1897

metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities

1664

for i in range(n):

1665

nMax += nLength

1666

phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))

1667

phaseCCF = numpy.mean(phaseCCF, axis = 2)

1668

phaseOffset[:,i] = phaseCCF.transpose()

1669

nMin = nMax

1670

# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)

1898

1671

1899

velEst = numpy.zeros((heightList.size,2))*numpy.nan

1672

#Remove Outliers

1900

azimuth1 = azimuth1*numpy.pi/180

1673

factor = 2

1674

wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))

1675

dw = numpy.std(wt,axis = 1)

1676

dw = dw.reshape((dw.size,1))

1677

ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))

1678

phaseOffset[ind] = numpy.nan

1679

phaseOffset = stats.nanmean(phaseOffset, axis=1)

1901

1680

1902

for i in range(heightList.size):

1681

return phaseOffset

1903

h = heightList[i]

1904

indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]

1905

metHeight = metArray1[indH,:]

1906

if metHeight.shape[0] >= 2:

1907

velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities

1908

iazim = metHeight[:,1].astype(int)

1909

azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths

1910

A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))

1911

A = numpy.asmatrix(A)

1912

A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()

1913

velHor = numpy.dot(A1,velAux)

1914

1915

velEst[i,:] = numpy.squeeze(velHor)

1916

return velEst

1917

1682

1918

def __getPhase~~Slope~~(self, ~~metArray~~, ~~heightList~~, ~~timeLis~~t):

1683

def __shiftPhase(self, data, phaseShift):

1919

meteorList = []

1684

#this will shift the phase of a complex number

1920

#utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2

1685

dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)

1921

#Putting back together the meteor matrix

1686

return dataShifted

1922

utctime = metArray[:,0]

1687

1923

uniqueTime = numpy.unique(utctime)

1688

def __estimatePhaseDifference(self, array, pairslist):

1924

1689

nChannel = array.shape[0]

1925

phaseDerThresh = 0.5

1690

nHeights = array.shape[2]

1926

ippSeconds = timeList[1] - timeList[0]

1691

numPairs = len(pairslist)

1927

sec = numpy.where(timeList>1)[0][0]

1692

# phaseCCF = numpy.zeros((nChannel, 5, nHeights))

1928

nPairs = metArray.shape[1] - 6

1693

phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))

1929

nHeights = len(heightList)

1930

1694

1931

for t in uniqueTime:

1695

#Correct phases

1932

metArray1 = metArray[utctime==t,:]

1696

derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]

1933

# phaseDerThresh = numpy.pi/4 #reducir Phase thresh

1697

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

1934

tmet = metArray1[:,1].astype(int)

1935

hmet = metArray1[:,2].astype(int)

1936

1937

metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))

1938

metPhase[:,:] = numpy.nan

1939

metPhase[:,hmet,tmet] = metArray1[:,6:].T

1940

1941

#Delete short trails

1942

metBool = ~numpy.isnan(metPhase[0,:,:])

1943

heightVect = numpy.sum(metBool, axis = 1)

1944

metBool[heightVect<sec,:] = False

1945

metPhase[:,heightVect<sec,:] = numpy.nan

1946

1947

#Derivative

1948

metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])

1949

phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))

1950

metPhase[phDerAux] = numpy.nan

1951

1952

#--------------------------METEOR DETECTION -----------------------------------------

1953

indMet = numpy.where(numpy.any(metBool,axis=1))[0]

1954

1955

for p in numpy.arange(nPairs):

1956

phase = metPhase[p,:,:]

1957

phDer = metDer[p,:,:]

1958

1959

for h in indMet:

1960

height = heightList[h]

1961

phase1 = phase[h,:] #82

1962

phDer1 = phDer[h,:]

1963

1964

phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap

1965

1966

indValid = numpy.where(~numpy.isnan(phase1))[0]

1967

initMet = indValid[0]

1968

endMet = 0

1969

1970

for i in range(len(indValid)-1):

1971

1972

#Time difference

1973

inow = indValid[i]

1974

inext = indValid[i+1]

1975

idiff = inext - inow

1976

#Phase difference

1977

phDiff = numpy.abs(phase1[inext] - phase1[inow])

1978

1979

if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor

1980

sizeTrail = inow - initMet + 1

1981

if sizeTrail>3*sec: #Too short meteors

1982

x = numpy.arange(initMet,inow+1)*ippSeconds

1983

y = phase1[initMet:inow+1]

1984

ynnan = ~numpy.isnan(y)

1985

x = x[ynnan]

1986

y = y[ynnan]

1987

slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)

1988

ylin = x*slope + intercept

1989

rsq = r_value**2

1990

if rsq > 0.5:

1991

vel = slope#*height*1000/(k*d)

1992

estAux = numpy.array([utctime,p,height, vel, rsq])

1993

meteorList.append(estAux)

1994

initMet = inext

1995

metArray2 = numpy.array(meteorList)

1996

1698

1997

return metArray2

1699

if indDer[0].shape[0] > 0:

1700

for i in range(indDer[0].shape[0]):

1701

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])

1702

phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi

1998

1703

1999

def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):

1704

# for j in range(numSides):

1705

# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])

1706

# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)

1707

#

1708

#Linear

1709

phaseInt = numpy.zeros((numPairs,1))

1710

angAllCCF = phaseCCF[:,[0,1,3,4],0]

1711

for j in range(numPairs):

1712

fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])

1713

phaseInt[j] = fit[1]

1714

#Phase Differences

1715

phaseDiff = phaseInt - phaseCCF[:,2,:]

1716

phaseArrival = phaseInt.reshape(phaseInt.size)

2000

1717

2001

azimuth1 = numpy.zeros(len(pairslist))

1718

#Dealias

2002

dist = numpy.zeros(len(pairslist))

1719

phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))

1720

# indAlias = numpy.where(phaseArrival > numpy.pi)

1721

# phaseArrival[indAlias] -= 2*numpy.pi

1722

# indAlias = numpy.where(phaseArrival < -numpy.pi)

1723

# phaseArrival[indAlias] += 2*numpy.pi

2003

1724

2004

for i in range(len(rx_location)):

1725

return phaseDiff, phaseArrival

2005

ch0 = pairslist[i][0]

1726

2006

ch1 = pairslist[i][1]

1727

def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):

1728

#this function will run the coherent detection used in Holdworth et al. 2004 and return the net power

1729

#find the phase shifts of each channel over 1 second intervals

1730

#only look at ranges below the beacon signal

1731

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

1732

numBlocks = int(volts.shape[1]/numProfPerBlock)

1733

numHeights = volts.shape[2]

1734

nChannel = volts.shape[0]

1735

voltsCohDet = volts.copy()

1736

1737

pairsarray = numpy.array(pairslist)

1738

indSides = pairsarray[:,1]

1739

# indSides = numpy.array(range(nChannel))

1740

# indSides = numpy.delete(indSides, indCenter)

1741

#

1742

# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)

1743

listBlocks = numpy.array_split(volts, numBlocks, 1)

1744

1745

startInd = 0

1746

endInd = 0

2007

1747

2008

diffX = rx_location[ch0][0] - rx_location[ch1][0]

1748

for i in range(numBlocks):

2009

diffY = rx_location[ch0][1] - rx_location[ch1][1]

1749

startInd = endInd

2010

azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi

1750

endInd = endInd + listBlocks[i].shape[1]

2011

dist[i] = numpy.sqrt(diffX**2 + diffY**2)

1751

1752

arrayBlock = listBlocks[i]

1753

# arrayBlockCenter = listCenter[i]

1754

1755

#Estimate the Phase Difference

1756

phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)

1757

#Phase Difference RMS

1758

arrayPhaseRMS = numpy.abs(phaseDiff)

1759

phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)

1760

indPhase = numpy.where(phaseRMSaux==4)

1761

#Shifting

1762

if indPhase[0].shape[0] > 0:

1763

for j in range(indSides.size):

1764

arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())

1765

voltsCohDet[:,startInd:endInd,:] = arrayBlock

1766

1767

return voltsCohDet

1768

1769

def __calculateCCF(self, volts, pairslist ,laglist):

2012

1770

2013

azimuth1 -= azimuth0

1771

nHeights = volts.shape[2]

2014

return azimuth1, dist

1772

nPoints = volts.shape[1]

1773

voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')

1774

1775

for i in range(len(pairslist)):

1776

volts1 = volts[pairslist[i][0]]

1777

volts2 = volts[pairslist[i][1]]

1778

1779

for t in range(len(laglist)):

1780

idxT = laglist[t]

1781

if idxT >= 0:

1782

vStacked = numpy.vstack((volts2[idxT:,:],

1783

numpy.zeros((idxT, nHeights),dtype='complex')))

1784

else:

1785

vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),

1786

volts2[:(nPoints + idxT),:]))

1787

voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)

2015

1788

2016

def techniqueNSM_DBS(self, **kwargs):

1789

vStacked = None

2017

metArray = kwargs['metArray']

1790

return voltsCCF

2018

heightList = kwargs['heightList']

2019

timeList = kwargs['timeList']

2020

zenithList = kwargs['zenithList']

2021

nChan = numpy.max(cmet) + 1

2022

nHeights = len(heightList)

2023

1791

2024

utctime = metArray[:,0]

1792

def __getNoise(self, power, timeSegment, timeInterval):

2025

cmet = metArray[:,1]

1793

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

2026

hmet = metArray1[:,3].astype(int)

1794

numBlocks = int(power.shape[0]/numProfPerBlock)

2027

h1met = heightList[hmet]*zenithList[cmet]

1795

numHeights = power.shape[1]

2028

vmet = metArray1[:,5]

1796

1797

listPower = numpy.array_split(power, numBlocks, 0)

1798

noise = numpy.zeros((power.shape[0], power.shape[1]))

1799

noise1 = numpy.zeros((power.shape[0], power.shape[1]))

2029

1800

2030

for i in range(nHeights - 1):

1801

startInd = 0

2031

hmin = heightList[i]

1802

endInd = 0

2032

hmax = heightList[i + 1]

1803

1804

for i in range(numBlocks): #split por canal

1805

startInd = endInd

1806

endInd = endInd + listPower[i].shape[0]

2033

1807

2034

vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]

1808

arrayBlock = listPower[i]

1809

noiseAux = numpy.mean(arrayBlock, 0)

1810

# noiseAux = numpy.median(noiseAux)

1811

# noiseAux = numpy.mean(arrayBlock)

1812

noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux

2035

1813

1814

noiseAux1 = numpy.mean(arrayBlock)

1815

noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1

1816

1817

return noise, noise1

1818

1819

def __findMeteors(self, power, thresh):

1820

nProf = power.shape[0]

1821

nHeights = power.shape[1]

1822

listMeteors = []

1823

1824

for i in range(nHeights):

1825

powerAux = power[:,i]

1826

threshAux = thresh[:,i]

2036

1827

1828

indUPthresh = numpy.where(powerAux > threshAux)[0]

1829

indDNthresh = numpy.where(powerAux <= threshAux)[0]

2037

1830

2038

return data_output

1831

j = 0

2039

1832

2040

def run(self, dataOut, technique, **kwargs):

1833

while (j < indUPthresh.size - 2):

1834

if (indUPthresh[j + 2] == indUPthresh[j] + 2):

1835

indDNAux = numpy.where(indDNthresh > indUPthresh[j])

1836

indDNthresh = indDNthresh[indDNAux]

1837

1838

if (indDNthresh.size > 0):

1839

indEnd = indDNthresh[0] - 1

1840

indInit = indUPthresh[j]

1841

1842

meteor = powerAux[indInit:indEnd + 1]

1843

indPeak = meteor.argmax() + indInit

1844

FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))

1845

1846

listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!

1847

j = numpy.where(indUPthresh == indEnd)[0] + 1

1848

else: j+=1

1849

else: j+=1

1850

1851

return listMeteors

1852

1853

def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):

2041

1854

2042

param = dataOut.data_param

1855

arrayMeteors = numpy.asarray(listMeteors)

2043

if dataOut.abscissaList != None:

1856

listMeteors1 = []

2044

absc = dataOut.abscissaList[:-1]

2045

noise = dataOut.noise

2046

heightList = dataOut.heightList

2047

SNR = dataOut.data_SNR

2048

1857

2049

if technique == 'DBS':

1858

while arrayMeteors.shape[0] > 0:

1859

FLAs = arrayMeteors[:,4]

1860

maxFLA = FLAs.argmax()

1861

listMeteors1.append(arrayMeteors[maxFLA,:])

2050

1862

2051

if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):

1863

MeteorInitTime = arrayMeteors[maxFLA,1]

2052

theta_x = numpy.array(kwargs['dirCosx'])

1864

MeteorEndTime = arrayMeteors[maxFLA,3]

2053

theta_y = numpy.array(kwargs['dirCosy'])

1865

MeteorHeight = arrayMeteors[maxFLA,0]

2054

else:

2055

elev = numpy.array(kwargs['elevation'])

2056

azim = numpy.array(kwargs['azimuth'])

2057

theta_x, theta_y = self.__calculateCosDir(elev, azim)

2058

azimuth = kwargs['correctAzimuth']

2059

if kwargs.has_key('horizontalOnly'):

2060

horizontalOnly = kwargs['horizontalOnly']

2061

else: horizontalOnly = False

2062

if kwargs.has_key('correctFactor'):

2063

correctFactor = kwargs['correctFactor']

2064

else: correctFactor = 1

2065

if kwargs.has_key('channelList'):

2066

channelList = kwargs['channelList']

2067

if len(channelList) == 2:

2068

horizontalOnly = True

2069

arrayChannel = numpy.array(channelList)

2070

param = param[arrayChannel,:,:]

2071

theta_x = theta_x[arrayChannel]

2072

theta_y = theta_y[arrayChannel]

2073

1866

2074

velRadial0 = param[:,1,:] #Radial velocity

1867

#Check neighborhood

2075

dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function

1868

maxHeightIndex = MeteorHeight + rangeLimit

2076

dataOut.utctimeInit = dataOut.utctime

1869

minHeightIndex = MeteorHeight - rangeLimit

2077

dataOut.outputInterval = dataOut.paramInterval

1870

minTimeIndex = MeteorInitTime - timeLimit

1871

maxTimeIndex = MeteorEndTime + timeLimit

2078

1872

2079

elif technique == 'SA':

1873

#Check Heights

1874

indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)

1875

indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)

1876

indBoth = numpy.where(numpy.logical_and(indTime,indHeight))

1877

1878

arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)

2080

1879

2081

#Parameters

1880

return listMeteors1

2082

position_x = kwargs['positionX']

1881

2083

position_y = kwargs['positionY']

1882

def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):

2084

azimuth = kwargs['azimuth']

1883

numHeights = volts.shape[2]

1884

nChannel = volts.shape[0]

1885

1886

thresholdPhase = thresh[0]

1887

thresholdNoise = thresh[1]

1888

thresholdDB = float(thresh[2])

1889

1890

thresholdDB1 = 10**(thresholdDB/10)

1891

pairsarray = numpy.array(pairslist)

1892

indSides = pairsarray[:,1]

1893

1894

pairslist1 = list(pairslist)

1895

pairslist1.append((0,1))

1896

pairslist1.append((3,4))

1897

1898

listMeteors1 = []

1899

listPowerSeries = []

1900

listVoltageSeries = []

1901

#volts has the war data

1902

1903

if frequency == 30e6:

1904

timeLag = 45*10**-3

1905

else:

1906

timeLag = 15*10**-3

1907

lag = numpy.ceil(timeLag/timeInterval)

1908

1909

for i in range(len(listMeteors)):

1910

1911

###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################

1912

meteorAux = numpy.zeros(16)

1913

1914

#Loading meteor Data (mHeight, mStart, mPeak, mEnd)

1915

mHeight = listMeteors[i][0]

1916

mStart = listMeteors[i][1]

1917

mPeak = listMeteors[i][2]

1918

mEnd = listMeteors[i][3]

1919

1920

#get the volt data between the start and end times of the meteor

1921

meteorVolts = volts[:,mStart:mEnd+1,mHeight]

1922

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

2085

1923

2086

if kwargs.has_key('crosspairsList'):

1924

#3.6. Phase Difference estimation

2087

pairs = kwargs['crosspairsList']

1925

phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)

2088

else:

2089

pairs = None

2090

2091

if kwargs.has_key('correctFactor'):

2092

correctFactor = kwargs['correctFactor']

2093

else:

2094

correctFactor = 1

2095

2096

tau = dataOut.data_param

2097

_lambda = dataOut.C/dataOut.frequency

2098

pairsList = dataOut.groupList

2099

nChannels = dataOut.nChannels

2100

2101

dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)

2102

dataOut.utctimeInit = dataOut.utctime

2103

dataOut.outputInterval = dataOut.timeInterval

2104

1926

2105

elif technique == 'Meteors':

1927

#3.7. Phase difference removal & meteor start, peak and end times reestimated

2106

dataOut.flagNoData = True

1928

#meteorVolts0.- all Channels, all Profiles

2107

self.__dataReady = False

1929

meteorVolts0 = volts[:,:,mHeight]

1930

meteorThresh = noise[:,mHeight]*thresholdNoise

1931

meteorNoise = noise[:,mHeight]

1932

meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting

1933

powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power

2108

1934

2109

if kwargs.has_key('nHours'):

1935

#Times reestimation

2110

nHours = kwargs['nHours']

1936

mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]

2111

else:

1937

if mStart1.size > 0:

2112

~~nHours~~ = 1

1938

mStart1 = mStart1[-1] + 1

2113

2114

if kwargs.has_key('meteorsPerBin'):

2115

meteorThresh = kwargs['meteorsPerBin']

2116

else:

2117

meteorThresh = 6

2118

1939

2119

if kwargs.has_key('hmin'):

1940

else:

2120

~~hmin~~ = k~~wargs~~[~~'hmin'~~]

1941

mStart1 = mPeak

2121

else: hmin = 70

2122

if kwargs.has_key('hmax'):

2123

hmax = kwargs['hmax']

2124

else: hmax = 110

2125

2126

dataOut.outputInterval = nHours*3600

2127

2128

if self.__isConfig == False:

2129

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

2130

#Get Initial LTC time

2131

self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

2132

self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

2133

2134

self.__isConfig = True

2135

1942

2136

if self.__buffer == None:

1943

mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1

2137

self.__buffer = dataOut.data_param

1944

mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]

2138

self.__firstdata = copy.copy(dataOut)

1945

if mEndDecayTime1.size == 0:

2139

1946

mEndDecayTime1 = powerNet0.size

2140

else:

1947

else:

2141

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

1948

mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1

1949

# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()

2142

1950

2143

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

1951

#meteorVolts1.- all Channels, from start to end

1952

meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]

1953

meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]

1954

if meteorVolts2.shape[1] == 0:

1955

meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]

1956

meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)

1957

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)

1958

##################### END PARAMETERS REESTIMATION #########################

2144

1959

2145

if self.__dataReady:

1960

##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################

2146

dataOut.utctimeInit = self.__initime

1961

# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis

1962

if meteorVolts2.shape[1] > 0:

1963

#Phase Difference re-estimation

1964

phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation

1965

# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)

1966

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])

1967

phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))

1968

meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting

2147

1969

2148

self.__initime += dataOut.outputInterval #to erase time offset

1970

#Phase Difference RMS

1971

phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))

1972

powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)

1973

#Data from Meteor

1974

mPeak1 = powerNet1.argmax() + mStart1

1975

mPeakPower1 = powerNet1.max()

1976

noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])

1977

mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux

1978

Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])

1979

Meteor1 = numpy.hstack((Meteor1,phaseDiffint))

1980

PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]

1981

#Vectorize

1982

meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]

1983

meteorAux[7:11] = phaseDiffint[0:4]

2149

1984

2150

dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)

1985

#Rejection Criterions

2151

dataOut.flagNoData = False

1986

if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation

2152

self.__buffer = None

1987

meteorAux[-1] = 17

1988

elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB

1989

meteorAux[-1] = 1

2153

1990

2154

elif technique == 'Meteors1':

2155

dataOut.flagNoData = True

2156

self.__dataReady = False

2157

2158

if kwargs.has_key('nMins'):

2159

nMins = kwargs['nMins']

2160

else: nMins = 20

2161

if kwargs.has_key('rx_location'):

2162

rx_location = kwargs['rx_location']

2163

else: rx_location = [(0,1),(1,1),(1,0)]

2164

if kwargs.has_key('azimuth'):

2165

azimuth = kwargs['azimuth']

2166

else: azimuth = 51

2167

if kwargs.has_key('dfactor'):

2168

dfactor = kwargs['dfactor']

2169

if kwargs.has_key('mode'):

2170

mode = kwargs['mode']

2171

else: mode = 'SA'

2172

1991

2173

#Borrar luego esto

1992

else:

2174

if dataOut.groupList == None:

1993

meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]

2175

dataOut.groupList = [(0,1),(0,2),(1,2)]

1994

meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis

2176

groupList = dataOut.groupList

1995

PowerSeries = 0

2177

C = ~~3e8~~

1996

2178

freq = 50e6

1997

listMeteors1.append(meteorAux)

2179

lamb = C/freq

1998

listPowerSeries.append(PowerSeries)

2180

k = 2*numpy.pi/lamb

1999

listVoltageSeries.append(meteorVolts1)

2000

2001

return listMeteors1, listPowerSeries, listVoltageSeries

2002

2003

def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):

2004

2005

threshError = 10

2006

#Depending if it is 30 or 50 MHz

2007

if frequency == 30e6:

2008

timeLag = 45*10**-3

2009

else:

2010

timeLag = 15*10**-3

2011

lag = numpy.ceil(timeLag/timeInterval)

2012

2013

listMeteors1 = []

2014

2015

for i in range(len(listMeteors)):

2016

meteorPower = listPower[i]

2017

meteorAux = listMeteors[i]

2181

2018

2182

timeList = dataOut.abscissaList

2019

if meteorAux[-1] == 0:

2183

heightList = dataOut.heightList

2020

2021

try:

2022

indmax = meteorPower.argmax()

2023

indlag = indmax + lag

2024

2025

y = meteorPower[indlag:]

2026

x = numpy.arange(0, y.size)*timeLag

2027

2028

#first guess

2029

a = y[0]

2030

tau = timeLag

2031

#exponential fit

2032

popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])

2033

y1 = self.__exponential_function(x, *popt)

2034

#error estimation

2035

error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))

2036

2037

decayTime = popt[1]

2038

riseTime = indmax*timeInterval

2039

meteorAux[11:13] = [decayTime, error]

2040

2041

#Table items 7, 8 and 11

2042

if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s

2043

meteorAux[-1] = 7

2044

elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time

2045

meteorAux[-1] = 8

2046

if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time

2047

meteorAux[-1] = 11

2048

2049

2050

except:

2051

meteorAux[-1] = 11

2052

2184

2053

2185

if self.__isConfig == False:

2054

listMeteors1.append(meteorAux)

2186

dataOut.outputInterval = nMins*60

2055

2187

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

2056

return listMeteors1

2188

#Get Initial LTC time

2189

initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

2190

minuteAux = initime.minute

2191

minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)

2192

self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

2193

2057

2194

self.__isConfig = True

2058

#Exponential Function

2059

2060

def __exponential_function(self, x, a, tau):

2061

y = a*numpy.exp(-x/tau)

2062

return y

2063

2064

def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):

2065

2066

pairslist1 = list(pairslist)

2067

pairslist1.append((0,1))

2068

pairslist1.append((3,4))

2069

numPairs = len(pairslist1)

2070

#Time Lag

2071

timeLag = 45*10**-3

2072

c = 3e8

2073

lag = numpy.ceil(timeLag/timeInterval)

2074

freq = 30e6

2075

2076

listMeteors1 = []

2077

2078

for i in range(len(listMeteors)):

2079

meteorAux = listMeteors[i]

2080

if meteorAux[-1] == 0:

2081

mStart = listMeteors[i][1]

2082

mPeak = listMeteors[i][2]

2083

mLag = mPeak - mStart + lag

2084

2085

#get the volt data between the start and end times of the meteor

2086

meteorVolts = listVolts[i]

2087

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

2088

2089

#Get CCF

2090

allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])

2091

2092

#Method 2

2093

slopes = numpy.zeros(numPairs)

2094

time = numpy.array([-2,-1,1,2])*timeInterval

2095

angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])

2096

2097

#Correct phases

2098

derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]

2099

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

2195

2100

2196

if self.__buffer == None:

2101

if indDer[0].shape[0] > 0:

2197

self.__buffer = dataOut.data_param

2102

for i in range(indDer[0].shape[0]):

2198

self.__firstdata = copy.copy(dataOut)

2103

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])

2104

angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi

2199

2105

2200

else:

2106

# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))

2201

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

2107

for j in range(numPairs):

2202

2108

fit = stats.linregress(time, angAllCCF[j,:])

2203

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

2109

slopes[j] = fit[0]

2204

2205

if self.__dataReady:

2206

dataOut.utctimeInit = self.__initime

2207

self.__initime += dataOut.outputInterval #to erase time offset

2208

2110

2209

~~metArray~~ = ~~self~~.~~__buff~~er

2111

#Remove Outlier

2210

if mode == 'SA':

2112

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

2211

dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)

2113

# slopes = numpy.delete(slopes,indOut)

2212

elif mode == 'DBS':

2114

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

2213

dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)

2115

# slopes = numpy.delete(slopes,indOut)

2214

dataOut.data_output = dataOut.data_output.T

2116

2215

dataOut.flagNoData = False

2117

radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)

2216

self.__buffer = None

2118

radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)

2217

2119

meteorAux[-2] = radialError

2218

return

2120

meteorAux[-3] = radialVelocity

2219

2121

2220

class EWDriftsEstimation(Operation):

2122

#Setting Error

2221

2123

#Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s

2222

def __init__(self):

2124

if numpy.abs(radialVelocity) > 200:

2223

Operation.__init__(self)

2125

meteorAux[-1] = 15

2126

#Number 12: Poor fit to CCF variation for estimation of radial drift velocity

2127

elif radialError > radialStdThresh:

2128

meteorAux[-1] = 12

2129

2130

listMeteors1.append(meteorAux)

2131

return listMeteors1

2224

2132

2225

def __correctValues(self, heiRang, phi, velRadial, SNR):

2133

def __setNewArrays(self, listMeteors, date, heiRang):

2226

listPhi = phi.tolist()

2227

maxid = listPhi.index(max(listPhi))

2228

minid = listPhi.index(min(listPhi))

2229

2134

2230

rango = range(len(phi))

2135

#New arrays

2231

# rango = numpy.delete(rango,maxid)

2136

arrayMeteors = numpy.array(listMeteors)

2137

arrayParameters = numpy.zeros((len(listMeteors), 13))

2232

2138

2233

heiRang1 = heiRang*math.cos(phi[maxid])

2139

#Date inclusion

2234

heiRangAux = heiRang*math.cos(phi[minid])

2140

# date = re.findall(r'\((.*?)\)', date)

2235

indOut = (heiRang1 < heiRangAux[0]).nonzero()

2141

# date = date[0].split(',')

2236

heiRang1 = numpy.delete(heiRang1,indOut)

2142

# date = map(int, date)

2143

#

2144

# if len(date)<6:

2145

# date.append(0)

2146

#

2147

# date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]

2148

# arrayDate = numpy.tile(date, (len(listMeteors), 1))

2149

arrayDate = numpy.tile(date, (len(listMeteors)))

2237

2150

2238

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

2151

#Meteor array

2239

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

2152

# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]

2153

# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))

2240

2154

2241

for i in rango:

2155

#Parameters Array

2242

x = heiRang*math.cos(phi[i])

2156

arrayParameters[:,0] = arrayDate #Date

2243

y1 = velRadial[i,:]

2157

arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range

2244

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

2158

arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error

2245

2159

arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases

2246

x1 = heiRang1

2160

arrayParameters[:,-1] = arrayMeteors[:,-1] #Error

2247

y11 = f1(x1)

2248

2249

y2 = SNR[i,:]

2250

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

2251

y21 = f2(x1)

2252

2253

velRadial1[i,:] = y11

2254

SNR1[i,:] = y21

2255

2256

return heiRang1, velRadial1, SNR1

2257

2161

2258

def run(self, dataOut, zenith, zenithCorrection):

2259

heiRang = dataOut.heightList

2260

velRadial = dataOut.data_param[:,3,:]

2261

SNR = dataOut.data_SNR

2262

2162

2263

zenith = numpy.array(zenith)

2163

return arrayParameters

2264

zenith -= zenithCorrection

2164

2265

zenith *= numpy.pi/180

2165

class CorrectMeteorPhases(Operation):

2266

2166

2267

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)

2167

def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):

2268

2168

2269

alp = zenith[0]

2169

arrayParameters = dataOut.data_param

2270

bet = ~~zenith~~[1]

2170

pairsList = []

2271

2171

pairx = (0,1)

2272

w_w = velRadial1[0,:]

2172

pairy = (2,3)

2273

w_e = velRadial1[1,:]

2173

pairsList.append(pairx)

2274

2174

pairsList.append(pairy)

2275

w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))

2175

jph = numpy.zeros(4)

2276

u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))

2277

2176

2278

winds = numpy.vstack((u,w))

2177

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

2178

# arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

2179

arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))

2279

2180

2280

dataOut.heightList = heiRang1

2181

meteorOps = MeteorOperations()

2281

dataOut.data_output = winds

2182

if channelPositions == None:

2282

dataOut.data_SNR = SNR1

2183

# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T

2184

channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella

2185

2186

pairslist0, distances = meteorOps.getPhasePairs(channelPositions)

2187

h = (hmin,hmax)

2188

2189

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

2283

2190

2284

dataOut.utctimeInit = dataOut.utctime

2191

dataOut.data_param = arrayParameters

2285

dataOut.outputInterval = dataOut.timeInterval

2286

return

2192

return

2287

2193

2288

class PhaseCalibration(Operation):

2194

class PhaseCalibration(Operation):

2289

2195

2290

__buffer = None

2196

__buffer = None

2709

# pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]

2615

# pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]

2710

pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]

2616

pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]

2711

2617

2712

return pairslist, distances No newline at end of file

2618

return pairslist, distances

2619

# def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):

2620

#

2621

# arrayAOA = numpy.zeros((phases.shape[0],3))

2622

# cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)

2623

#

2624

# arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)

2625

# cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)

2626

# arrayAOA[:,2] = cosDirError

2627

#

2628

# azimuthAngle = arrayAOA[:,0]

2629

# zenithAngle = arrayAOA[:,1]

2630

#

2631

# #Setting Error

2632

# #Number 3: AOA not fesible

2633

# indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]

2634

# error[indInvalid] = 3

2635

# #Number 4: Large difference in AOAs obtained from different antenna baselines

2636

# indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]

2637

# error[indInvalid] = 4

2638

# return arrayAOA, error

2639

#

2640

# def __getDirectionCosines(self, arrayPhase, pairsList):

2641

#

2642

# #Initializing some variables

2643

# ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi

2644

# ang_aux = ang_aux.reshape(1,ang_aux.size)

2645

#

2646

# cosdir = numpy.zeros((arrayPhase.shape[0],2))

2647

# cosdir0 = numpy.zeros((arrayPhase.shape[0],2))

2648

#

2649

#

2650

# for i in range(2):

2651

# #First Estimation

2652

# phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]

2653

# #Dealias

2654

# indcsi = numpy.where(phi0_aux > numpy.pi)

2655

# phi0_aux[indcsi] -= 2*numpy.pi

2656

# indcsi = numpy.where(phi0_aux < -numpy.pi)

2657

# phi0_aux[indcsi] += 2*numpy.pi

2658

# #Direction Cosine 0

2659

# cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)

2660

#

2661

# #Most-Accurate Second Estimation

2662

# phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]

2663

# phi1_aux = phi1_aux.reshape(phi1_aux.size,1)

2664

# #Direction Cosine 1

2665

# cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)

2666

#

2667

# #Searching the correct Direction Cosine

2668

# cosdir0_aux = cosdir0[:,i]

2669

# cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)

2670

# #Minimum Distance

2671

# cosDiff = (cosdir1 - cosdir0_aux)**2

2672

# indcos = cosDiff.argmin(axis = 1)

2673

# #Saving Value obtained

2674

# cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]

2675

#

2676

# return cosdir0, cosdir

2677

#

2678

# def __calculateAOA(self, cosdir, azimuth):

2679

# cosdirX = cosdir[:,0]

2680

# cosdirY = cosdir[:,1]

2681

#

2682

# zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi

2683

# azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east

2684

# angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()

2685

#

2686

# return angles

2687

#

2688

# def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):

2689

#

2690

# Ramb = 375 #Ramb = c/(2*PRF)

2691

# Re = 6371 #Earth Radius

2692

# heights = numpy.zeros(Ranges.shape)

2693

#

2694

# R_aux = numpy.array([0,1,2])*Ramb

2695

# R_aux = R_aux.reshape(1,R_aux.size)

2696

#

2697

# Ranges = Ranges.reshape(Ranges.size,1)

2698

#

2699

# Ri = Ranges + R_aux

2700

# hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re

2701

#

2702

# #Check if there is a height between 70 and 110 km

2703

# h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)

2704

# ind_h = numpy.where(h_bool == 1)[0]

2705

#

2706

# hCorr = hi[ind_h, :]

2707

# ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))

2708

#

2709

# hCorr = hi[ind_hCorr]

2710

# heights[ind_h] = hCorr

2711

#

2712

# #Setting Error

2713

# #Number 13: Height unresolvable echo: not valid height within 70 to 110 km

2714

# #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km

2715

#

2716

# indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]

2717

# error[indInvalid2] = 14

2718

# indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]

2719

# error[indInvalid1] = 13

2720

#

2721

# return heights, error

2722

No newline at end of file

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages