'LoranView 'DF6NM, 2004-2020 '1q20 200618: gri 6731 jumps to track frequency offset '1q19 200511: remove 8000 from pps aquisition '1q18 180620: bugfix gri5990step txdelayaqA(nst) '1q17 171027/28: remove nextpulsestart# drift due to rounding errors '1q16a 170218: gri 5990 jumps also to txdelayA (ms) '1q15 170212: gri 5990 daily jumps '1q14 150826: long path with l flag (e.g. 8970m l) '1q13 150330: nb intensity correction '1q12a 150330: notch100khz buffer averages after nb '1q12 150329: nb apow!=0 after loran blanking, '1q11 150325: nbspread, nbthre! reset on resync, nb=2 '1q10 150323: noise blanker '1q9 150318: notch100khz bugfix ' (todo: all-in-GRI sync, pulsetime rounding errors) '1q8 141017: TAI guess 1pps nearest to PC clock, select bmp and table filenames '1q7 141014: sync to nearest 2GRI (startoffs# < t2gri0#/2) '1q6 141008: pctime! corrected for 1 channel '1q5 141007: bmpcmdcount, statuslines with linefeeds '1q4 141004: bmp fileaccess conflict '1q3 140929: timer midnight bugfix, crossrate cancellation xrc7950 '1q2 140929: subsecond timer '1q1 140929: some code cleanup, leapsec 120630 fixed, sleep (needs FreeBasic) '1p11 140928: srcorrpfac! '1p10 140923: 15000 km Vincenty atan2, taisecs! '1p9 140325: FreeBasic: fbc /lang qb -ex lv_1p9.bas '1p8 120906: 8000->6000 crossrate cancellation '1p7 120903: edge search from left to right, lastLOus '1p6 120716: leapsec 120630, clockoffset '1p5 120714: lo freq output '1p4 120522: lo drift '1p3 110409: lo pll cleanup '1p2 110408: lo pll with freq-to-voltage-to-varicap '1p1 110408: lo pll using sndoutpt and lo/8 coherent audio '1o5 110407: lo phase time corrected for sr drift '1o4 111219: fixed 24h pc-sr time '1o3 111218: blankingdistance, bugfixes for sync with slaves, bmp in use '1o2 111214: eliminated cormr cormi buffers, added clockerrorlimit, clockoffset '1o1 111213: GPS sync, bugfixes: corend, nmsssm=4 (1 master + 3 slaves) '1n9 111213: accurate PC clock time for sync '1n8 100316: dat=other filename for SndInput '1n7 100209: allow masters KLMN slaves OPQRSTUVWXYZ '1n6 100203: leapsecs 2008/09, 100 kHz notch '1n5 100202: fixed out of heap space in interpreter mode '1n4 090209: filescale '1n3 090209: fileblocklen (256 required for eurofix) '1n2 090129: sleeptime, searchrange '1n1 080908: vincenty ellipsoid, nedm=20, dds lo, hh:mm:ss 080913: bug fixes '1m6 060318: fast LO > sr/32, individual pulse phase for 137 kHz monitoring, srcorfac! '1m5 060314: west-to-east sorted display '1m4 060312: debugged fast LO, CRLF for aquisition time, name first in stations.txt, ASF > 1 '1m3 060130: selectable ref station '1m2 060127: ' variable blanking interval, slightly faster averaging loop, ' newline printed , fixed nominal sr for looffset '1m1 060107: faster looffset mixer '1l9 '1e 040120: ... '1 040117: first averaging '$ lang: "qb" DEFLNG A-Z PRINT PRINT " >>> LoranView (1q20) by DF6NM - press q to quit <<<" PRINT 'initialize: inifile$ = "loranvw.ini" cmd$ = COMMAND$ IF LEN(cmd$) > 0 THEN inifile$ = cmd$ 'DF6NM, nuernberg langwasser, jn59nj mylat# = 49.415762#: mylon# = 11.135018# wavflag = 0 file1$ = "audio.dat" file2$ = "lorvw.bmp" tablefile$ = "table.txt" date1$ = DATE$: time1$ = TIME$ utoffsh = 0 sr# = 22050 'sr# = 22042.453# 'Amilo 'sr# = 22050.7# '233MHz Notebook pi# = 4# * ATN(1#): pi! = pi# asf! = 1! 'sea water groundwave asf! = 1.008 'average for nighttime skywave 'c! = 299792458# / 1.0009 / asf! 'incl. primary and secondary factors for seawater shortintegrationtime! = 30 'seconds longintegrationfactor = 10 phasecorrflag = 1 srcorrfac! = 1 srcorrpfac! = 1 displaysamples = 24 centralsample = 6 firstblankedsample% = centralsample - 3 lastblankedsample% = centralsample + 4 'was + 5 - 1 resynclimit! = 0 '20 blankedstns = 0 hidden = 0 singlepulses = 0 'normal 18/16 pulse scheme colourflag = 1 fullscale! = 200 'LSB gamma! = 1 dbrange! = 0 linlogflag = 0 normalizedistance! = 0 normalizepower! = 0 hipass! = .1 hipassshift! = 1 edgelevel! = .4 referencephase! = 0 nstref = 0 blankingoffsetcorr! = 1 sleeptime! = 0 searchrange! = 3700 blocklen = 0 'default filescale! = 1! bytespersample% = 2 channels% = 2 bmperrorhandling%=0 'bmpcmd$="copy >nul lorvw.bmp lorvw2.bmp /y" bmpcmd$="" bmpcmdcount = 1 'call bmpcmd after every new bmp line lofreq# = 0 lodrift! = 0 locycles! = 0 hoursnow = -1 bandcenter# = 0 notch100kHzflag% = 1 xrcdelay7950! = 1.06 'ms n*53/50kHz cancels 6731 clockerrorlimit! = 0 'max clock error (0 = inactive) clockoffset! = 0 'PC clock error seconds + leap seconds after 2014 timercalavfac!=.01 timercalcnt=0 syncverbose% = 0 blankingdistance! = 0'km, 0 = inaktive pllflag% = 0 pllcmd$ = "SndOutp1 /testsig=" pllcenter# = 900 pllrange! = 850 pllfactort! = 10 pllfactorf! = 20 * 1E+08 plldither! = 0 pllavf! = .1 nbflag%=0 nbinit! = 100^2 'initial threshold nbfac!=2 ' threshold / average power nbupdate!=.1 'nb threshold adaption per shortintegration nbspread%=2 'high threshold at neighbouring pixels nbcorr!=1 'intensity correction factor gri5990step#=0 gri6731step#=0 remove8000frompps%=0 'readini: nstaqm = 10: nstm = 200: nedm = 20: nfirm = 20 'allow more chain offsets DIM griaq1A(nstaqm - 1), mxyzaq1$(nstaqm - 1), weightaq1!(nstaqm - 1) DIM gri1A(nstm - 1), mxyz1$(nstm - 1), stnflags$(nstm-1) DIM gried1(nedm - 1), ed1!(nedm - 1) DIM afirr!(nfirm - 1), afiri!(nfirm - 1) nstaq = 0: nst = 0: ned = 0: nfir = 0 OPEN inifile$ FOR INPUT AS 1 WHILE NOT EOF(1) LINE INPUT #1, in$ IF INSTR(in$, ";") > 0 THEN in$ = LEFT$(in$, INSTR(in$, ";") - 1) IF INSTR(in$, "'") > 0 THEN in$ = LEFT$(in$, INSTR(in$, "'") - 1) in$ = LCASE$(LTRIM$(RTRIM$(in$))) IF LEN(in$) > 0 THEN IF LEFT$(in$, 1) = "[" THEN topic$ = in$ topic$ = MID$(in$, 2, LEN(in$) - 2) ELSE equal = INSTR(in$, "=") IF equal > 0 THEN name$ = RTRIM$(LEFT$(in$, equal - 1)) val$ = LTRIM$((RIGHT$(in$, LEN(in$) - equal))) ELSE name$ = in$: val$ = "" END IF GOSUB evaluateiniline END IF END IF WEND CLOSE 1 nstaqm = nstaq nstm = nst nedm = ned nfirm = nfir DIM griaqA(nstaqm - 1), mxyzaq$(nstaqm - 1), weightaq!(nstaqm - 1) DIM deldraq!(nstaqm - 1) DIM griA(nstm - 1), mxyzA$(nstm - 1) DIM edaq!(nstaqm - 1), ed!(nstm - 1) DIM blankingflag%(nstm - 1) FOR n = 0 TO nstaqm - 1 griaqA(n) = griaq1A(n) mxyzaq$(n) = mxyzaq1$(n) weightaq!(n) = weightaq1!(n) FOR ned = 0 TO nedm - 1 IF griaqA(n) = gried1(ned) THEN edaq!(n) = ed1!(ned) NEXT ned NEXT n numblanked% = 0 FOR n = 0 TO nstm - 1 griA(n) = ABS(gri1A(n)) mxyzA$(n) = mxyz1$(n) IF (gri1A(n) < 0) AND numblanked% < blankedstns THEN blankingflag%(n) = 1 numblanked% = numblanked% + 1 END IF FOR ned = 0 TO nedm - 1 IF griA(n) = gried1(ned) THEN ed!(n) = ed1!(ned) NEXT ned NEXT n ERASE griaq1A, mxyzaq1$, weightaq1!, gri1A, mxyz1$, gried1, ed1! nstdm = nstm - hidden 'readtable: DIM location$(nstm - 1), pwrA!(nstm - 1), azimut!(nstm - 1) DIM txdelayA!(nstm - 1), dist!(nstm - 1) DIM txdelayaqA!(nstaqm - 1), distaq!(nstaqm - 1) ntab = 0 OPEN tablefile$ FOR INPUT AS 1 WHILE NOT EOF(1) INPUT #1, gri, gri2, mxyz$, name$ INPUT #1, latdeg, latmin, latsec!, ns$, londeg, lonmin, lonsec!, ew$ INPUT #1, txdelay!, codedelay!, pwr! lat# = ABS(latdeg) + ABS(latmin) / 60# + ABS(latsec!) / 3600# IF UCASE$(ns$) = "S" THEN lat# = -lat# lon# = ABS(londeg) + ABS(lonmin) / 60# + ABS(lonsec!) / 3600# IF UCASE$(ew$) = "W" THEN lon# = -lon# IF gri > 0 THEN mxyz$ = UCASE$(mxyz$) FOR nst = 0 TO nstm - 1 IF gri = griA(nst) AND mxyz$ = mxyzA$(nst) THEN 'GOSUB calculatedistance GOSUB calculatevincenty if instr(stnflags$(nst),"l")>0 then dist1#=dist# lat#=-lat#: lon#=lon#-180: if lon#<180 then lon#=lon#+360 GOSUB calculatevincenty dist#=dist#*2+dist1# lat#=-lat#: lon#=lon#+180: if lon#>180 then lon#=lon#-360 end if alf! = bet# '? location$(nst) = name$ txdelayA!(nst) = txdelay! / 1000! + ed!(nst) 'us -> ms pwrA!(nst) = pwr!: dist!(nst) = dist#: azimut!(nst) = alf! END IF NEXT nst FOR nst = 0 TO nstaqm - 1 IF gri = griaqA(nst) AND mxyz$ = mxyzaq$(nst) THEN 'GOSUB calculatedistance GOSUB calculatevincenty txdelayaqA!(nst) = txdelay! / 1000! + edaq!(nst) distaq!(nst) = dist# END IF NEXT nst ntab = ntab + 1 END IF WEND CLOSE 1 FOR nst = 0 TO nstm - 1 IF dist!(nst) < blankingdistance! AND numblanked% < blankedstns AND blankingflag%(nst) = 0 THEN blankingflag%(nst) = 1: numblanked% = numblanked% + 1 END IF NEXT nst FOR nst = 0 TO nstm - 1 PRINT griA(nst); mxyzA$(nst); IF dist!(nst) > 0 THEN azimut = azimut!(nst) PRINT " "; location$(nst); txdelayA!(nst); "ms"; PRINT pwrA!(nst); "kW"; dist!(nst); "km"; '; azimut; "deg"; ELSE PRINT " not found in table"; 'dist!(nst) = normalizedistance!: pwrA!(nst) = normalizepower! END IF IF nst >= nstdm THEN PRINT " (hidden)"; IF blankingflag%(nst) THEN PRINT " (blanked)"; PRINT stnflags$(nst) NEXT nst IF printstationlistflag% THEN OPEN "stations.txt" FOR OUTPUT AS 1 FOR nst = 0 TO nstdm - 1 PRINT #1, " "; location$(nst) PRINT #1, "_"; LTRIM$(RTRIM$(STR$(griA(nst)))); mxyzA$(nst); IF dist!(nst) > 0 THEN PRINT #1, RTRIM$(STR$(INT(dist!(nst) + .5))); "km" ELSE PRINT #1, " ?" END IF NEXT nst CLOSE 1 END IF PRINT numblanked%; "blanked, "; PRINT "sync to"; FOR nst = 0 TO nstaqm - 1 PRINT griaqA(nst); mxyzaq$(nst); NEXT nst DIM nst%(nstm - 1) n = 0 FOR nst = 0 TO nstm - 1 IF blankingflag%(nst) <> 0 THEN nst%(n) = nst: n = n + 1 NEXT nst FOR nst = 0 TO nstm - 1 IF blankingflag%(nst) = 0 THEN nst%(n) = nst: n = n + 1 NEXT nst 'start: inblocklen = 1024& * INT(sr# / 5512.5 + .5) 'inblocklen = 256& 'for Eurofix IF fileblocklen > 0 THEN inblocklen = fileblocklen bufoffs = 256& * (displaysamples \ 256& + 1) buflen = inblocklen + bufoffs inblockbytes = inblocklen * bytespersample% * channels% DIM ar%(buflen - 1), ai%(buflen - 1) DIM avr(nstm * displaysamples - 1), avi(nstm * displaysamples - 1) DIM avrA!(nstm * displaysamples - 1), aviA!(nstm * displaysamples - 1) DIM deldr!(nstm - 1), dtexp2#(nstm - 1) DIM avcntA(nstm - 1), serrorA!(nstm - 1), serrorl!(nstm - 1) DIM avcntsignA(nstm - 1) DIM nextpulsestartA#(nstm - 1), nextpulseingroup#(nstm - 1) DIM firstpulsestartA#(nstm - 1) '******** 1q17 DIM npulsemA(nstm - 1), npulseA(nstm - 1) DIM pulsespacing#(nstm - 1, 17), pulsesignA%(nstm - 1, 17) DIM amagA!(displaysamples - 1) DIM blankingfactorA!(nstm - 1) nmsssm = 4 'try to get master and up to 3 secondaries for sync correlation DIM aqbufstartA(nstaqm - 1), aqbufendA(nstaqm - 1), aqbufptrA(nstaqm - 1) DIM bufptr(nstaqm - 1) DIM startaq#(nstaqm - 1), avcntaq(nstaqm - 1) DIM amplA!(nstaqm, nmsssm - 1), delayA!(nstaqm - 1, nmsssm - 1) DIM dtcalcaq!(nstaqm - 1), dtmeasaq!(nstaqm - 1) DIM pulseoffscor(nstaqm - 1, 1, 17), pulsesigncor%(nstaqm - 1, 1, 17) DIM apow!(buflen-1), apow2!(nstm * displaysamples - 1), nbthr!(nstm * displaysamples - 1) DIM nbcnt!(nstm * displaysamples - 1) 'initialize pulsegroups sr0# = sr# searchpattern$ = "110010101" + "100111110" + "11111001" + "10101100" c! = 299792458# / 1.0009 / asf! FOR nst = 0 TO nstm - 1 IF INSTR("OPQRSTUVWXYZ", mxyzA$(nst)) > 0 AND singlepulses = 0 THEN npulsem = 16 FOR np = 0 TO 14 pulsespacing#(nst, np) = .001# NEXT np pulsespacing#(nst, 7) = .00001# * (griA(nst) - 700) pulsespacing#(nst, 15) = .00001# * (griA(nst) - 700) FOR np = 0 TO npulsem - 1 pulsesignA%(nst, np) = 2 * VAL(MID$(searchpattern$, 19 + np, 1)) - 1 NEXT np ELSE IF INSTR("KLMN", mxyzA$(nst)) > 0 AND singlepulses = 0 THEN npulsem = 18 FOR np = 0 TO 16 pulsespacing#(nst, np) = .001# NEXT np pulsespacing#(nst, 7) = .002# pulsespacing#(nst, 8) = .00001# * (griA(nst) - 900) pulsespacing#(nst, 16) = .002# pulsespacing#(nst, 17) = .00001# * (griA(nst) - 900) FOR np = 0 TO npulsem - 1 pulsesignA%(nst, np) = 2 * VAL(MID$(searchpattern$, 1 + np, 1)) - 1 NEXT np ELSE npulsem = 1 IF mxyzA$(nst) = "B" THEN npulsem = 2 IF singlepulses THEN npulsem = singlepulses pulsesign% = 1 FOR np = 0 TO npulsem - 1 pulsespacing#(nst, np) = .00002# * griA(nst) / npulsem pulsesignA%(nst, np) = pulsesign% pulsesign% = -pulsesign% NEXT np END IF END IF npulsemA(nst) = npulsem NEXT nst IF firstblankedsample% < 0 THEN firstblankedsample% = 0 IF lastblankedsample% >= displaysamples THEN lastblankedsample% = displaysamples - 1 blankedsamples% = lastblankedsample% - firstblankedsample% + 1 blankingfactor! = 1! FOR nsta% = 0 TO nstm - 1 nst = nst%(nsta%) blankingfactorA!(nst) = blankingfactor! IF blankingflag%(nst) THEN blankingfactor! = blankingfactor! * (1! - npulsemA(nst) * blankedsamples% / sr# / griA(nst) * 50000!) END IF NEXT nsta% DIM looffsetaq#(nstaqm - 1) bandcenterkhz = INT(bandcenter# / 1000 + .5) FOR naq = 0 TO nstaqm - 1 nearestgriharmonic = bandcenterkhz / 100 * griaqA(naq) looffsetaq#(naq) = nearestgriharmonic * 100000 / griaqA(naq) - bandcenter# NEXT naq bufstart# = -bufoffs / sr#: bufend# = 0 actualsamples = 0 srnominal! = sr# wr! = filescale!: wi! = 0: phcnt% = 0 if wavflag = 0 then gosub calibratetimer GOSUB openfile GOSUB evaluatedateandtime GOSUB synchronize GOSUB writebmpheader laq: resynccnt = 0 validcnt = 0 for n = 0 to nstm * displaysamples - 1 nbthr!(n) = nbinit! next n xrccnt6000 = 0: xrccnt7950 = 0 dd79508000! = 2!*50!/50000! 'vs 8000 dd79506731! = 2!*53!/50000! 'vs 6731 FOR nst = 0 TO nstm - 1 t2gri# = griA(nst) / 50000# deldr!(nst) = txdelayA!(nst) / 1000! + dist!(nst) * 1000! / c! '======= IF griA(nst) = 6000 THEN IF xrccnt6000 = 1 THEN deldr!(nst) = deldr!(nst) - .04 IF xrccnt6000 = 2 THEN deldr!(nst) = deldr!(nst) + .04 xrccnt6000 = xrccnt6000 + 1 END IF IF griA(nst) = 7950 THEN if xrccnt7950>=3 then xrccnt2 = (xrccnt7950-3) AND 1 IF xrccnt2 = 0 THEN deldr!(nst) = deldr!(nst) - xrcdelay7950!*.001 IF xrccnt2 = 1 THEN deldr!(nst) = deldr!(nst) + xrcdelay7950!*.001 '?:?xrccnt7950;xrccnt2;nst;deldr!(nst) end if xrccnt7950=xrccnt7950+1 END IF '======== expdelgris# = (deldr!(nst) - starttime#) / t2gri# 'time at rx dtexp2#(nst) = (expdelgris# - INT(expdelgris#)) * t2gri# npstart# = dtexp2#(nst) - (bufstart# + bufoffs / sr#) npstart# = npstart# - t2gri# * INT(npstart# / t2gri#) nextpulsestartA#(nst) = npstart# + bufstart# + (bufoffs - centralsample) / sr# firstpulsestartA#(nst)=nextpulsestartA#(nst) '**** 1q17 npulseA(nst) = 0 avcntA(nst) = 0 serrorA!(nst) = 0 avcntsignA(nst) = 0 NEXT nst '====== xrcoffs=0 if xrccnt6000 =3 then xrcoffs = -2 'if xrccnt7950 = 15 then xrcoffs = xrcoffs - 12 if xrccnt7950 = 9 then xrcoffs = xrcoffs - 6 '?xrcoffs '====== phaserefr! = COS(phasereference! * pi! / 180!) phaserefi! = SIN(phasereference! * pi! / 180!) phasefacr! = phaserefr!: phasefaci! = phaserefi! IF bufend# > 0 THEN integrationstart! = bufend# GOTO l1 'buf already contains data END IF WHILE bytesleft > 0 GOSUB getblockofdata l1: FOR nsta% = 0 TO nstm - 1 nst = nst%(nsta%) nbflagb%=nbflag% if nbflag%=2 then nbflagb%= 1%-blankingflag%(nst) nextpulsestart# = nextpulsestartA#(nst) npulse = npulseA(nst) avcnt = avcntA(nst) avcntsign = avcntsignA(nst) serror! = serrorA!(nst) ind2start% = nst * displaysamples ind2end% = ind2start% + displaysamples - 1 'now add pulses until before end of buffer WHILE nextpulsestart# + displaysamples / sr# < bufend# ind1# = (nextpulsestart# - bufstart#) * sr# ind1% = ind1# '- deltast! 'for jitter testing serror! = serror! + ind1# - ind1% IF bandcenter# THEN pulsephase# = bandcenter# * nextpulsestart# pulsephase! = 2 * pi# * (pulsephase# - INT(pulsephase#)) cosp! = COS(pulsephase!) * pulsesignA%(nst, npulse) sinp! = SIN(pulsephase!) * pulsesignA%(nst, npulse) FOR ind2% = ind2start% TO ind2end% ar! = ar%(ind1%): ai! = ai%(ind1%) ar = ar! * cosp! - ai! * sinp! ai = ar! * sinp! + ai! * cosp! avr(ind2%) = avr(ind2%) + ar avi(ind2%) = avi(ind2%) + ai ind1% = ind1% + 1 NEXT ind2% ELSE if nbflagb% then pulsesign% = pulsesignA%(nst, npulse) FOR ind2% = ind2start% TO ind2end% apow2!(ind2%)=apow2!(ind2%)+apow!(ind1%) if apow!(ind1%) < nbthr!(ind2%) then if pulsesign% > 0 THEN avr(ind2%) = avr(ind2%) + ar%(ind1%) avi(ind2%) = avi(ind2%) + ai%(ind1%) else avr(ind2%) = avr(ind2%) - ar%(ind1%) avi(ind2%) = avi(ind2%) - ai%(ind1%) end if else nbcnt!(ind2%)=nbcnt!(ind2%)+1! end if ind1% = ind1% + 1 NEXT ind2% else IF pulsesignA%(nst, npulse) > 0 THEN FOR ind2% = ind2start% TO ind2end% avr(ind2%) = avr(ind2%) + ar%(ind1%) avi(ind2%) = avi(ind2%) + ai%(ind1%) ind1% = ind1% + 1 NEXT ind2% ELSE FOR ind2% = ind2start% TO ind2end% avr(ind2%) = avr(ind2%) - ar%(ind1%) avi(ind2%) = avi(ind2%) - ai%(ind1%) ind1% = ind1% + 1 NEXT ind2% END IF end if END IF IF blankingflag%(nst) THEN 'blanking -3...+5 ind1% = ind1% - displaysamples + firstblankedsample% FOR nav% = ind1% TO ind1% + blankedsamples% - 1 ar%(nav%) = 0: ai%(nav%) = 0 apow!(nav%) = 0! NEXT nav% END IF avcntsign = avcntsign + pulsesignA%(nst, npulse) nextpulsestart# = nextpulsestart# + pulsespacing#(nst, npulse) '-1e-8 '**** 1q17 test npulse = npulse + 1 IF npulse = npulsemA(nst) THEN '***** 1q17 realign nextpulsestart# npulse = 0 ngris#=int((nextpulsestart#-firstpulsestartA#(nst))/(griA(nst)/100000#)+.5) nextpulsestart#=firstpulsestartA#(nst)+griA(nst)/100000#*ngris# END IF avcnt = avcnt + 1 WEND nextpulsestartA#(nst) = nextpulsestart# npulseA(nst) = npulse avcntA(nst) = avcnt serrorA!(nst) = serror! avcntsignA(nst) = avcntsign NEXT nsta% IF bufend# > integrationstart! + shortintegrationtime! THEN integrationstart! = integrationstart! + shortintegrationtime! if nbflag% then FOR nst = 0 TO nstm - 1 avcnt! = avcntA(nst) IF avcnt! > 0 THEN n%=nst*displaysamples for n1%= 0 to displaysamples-1 apowm!=apow2!(n%) for nspr%=1 to nbspread% if n1%-nspr%>=0 then if apow2!(n%-nspr%)>apowm! then apowm!=apow2!(n%-nspr%) end if if n1%+nspr%apowm! then apowm!=apow2!(n%+nspr%) end if next nspr% nbthr!(n%) = nbthr!(n%)*(1!-nbupdate!) + nbfac!*nbupdate!*apowm!/avcnt! n%=n%+1 next n1% for n%= nst*displaysamples to (nst+1)*displaysamples-1 apow2!(n%)=0! next n% end if next nst end if GOSUB printtime GOSUB measureslope longintcnt = longintcnt + 1 IF longintcnt >= longintegrationfactor THEN 'deltast! = deltast! + .1: IF deltast! > .3 THEN deltast! = -.3 'jitter testing longintcnt = 0 GOSUB writebmpline FOR n = 0 TO nstm * displaysamples - 1 avrA!(n) = 0: aviA!(n) = 0 NEXT n FOR nst = 0 TO nstm - 1 serrorl!(nst) = 0 NEXT nst END IF END IF IF resynccnt >= 3 THEN resynccnt = 0 validcnt = 0 IF wavflag THEN assumedstart# = assumedstart# + bufend# ELSE date1$ = DATE$: time1$ = TIME$: GOSUB evaluatedateandtime END IF sr# = sr0# bufstart# = -bufoffs / sr#: bufend# = 0: actualsamples = 0 GOSUB synchronize GOTO laq END IF WEND PRINT CLOSE 1 'GOSUB printstatistics END '============================ evaluateiniline: SELECT CASE topic$ CASE "location" SELECT CASE name$ CASE "lat": mylat# = VAL(val$) CASE "lon": mylon# = VAL(val$) END SELECT CASE "file" SELECT CASE name$ CASE "bmp": file2$ = val$ 'default lorvw.bmp CASE "table": tablefile$ = val$ 'table.txt CASE "wav": file1$ = val$: wavflag = 1 CASE "dat": file1$ = val$: wavflag = 0 'audio.dat CASE "date": date1$ = val$ CASE "time": time1$ = val$ CASE "bits": IF VAL(val$) <= 8 THEN bytespersample% = 1 ELSE bytespersample% = 2 CASE "stereo": channels% = 2 CASE "mono": channels% = 1 CASE "biasreal": offsetr! = VAL(val$) CASE "biasimag": offseti! = VAL(val$) CASE "fileblocklen": fileblocklen = VAL(val$) CASE "filescale": filescale! = VAL(val$) case "bmpcmd": bmpcmd$ = val$ case "bmpcmdcount": bmpcmdcount = val(val$) '-n = every nth new hour, 0 = never, +n = every n lines case "bmperrorhandling": bmperrorhandling% = val(val$) '0,1,2 END SELECT CASE "frequency" SELECT CASE name$ CASE "sr": sr# = VAL(val$) CASE "looffset": lofreq# = VAL(val$) CASE "lodrift": lodrift! = VAL(val$) CASE "srcorr": srcorrfac! = VAL(val$): srcorrpfac! = SQR(ABS(srcorrfac!)) CASE "phasecorr": phasecorrflag = VAL(val$) CASE "jittercorr": hipassshift! = VAL(val$) CASE "bandcenter": bandcenter# = VAL(val$) CASE "srlofactor": srlofactor! = VAL(val$) CASE "srlooffset": srlooffset! = VAL(val$) END SELECT CASE "pll" SELECT CASE name$ CASE "pllcommand": pllcmd$ = val$: pllflag% = 1 CASE "pllcenter": pllcenter# = VAL(val$) CASE "plllow": plllow! = VAL(val$) CASE "pllhigh": pllhigh! = VAL(val$) CASE "pllfactort": pllfactort! = VAL(val$) CASE "pllfactorf": pllfactorf! = VAL(val$) CASE "plldither": plldither! = VAL(val$) CASE "pllaveragef": pllavf! = VAL(val$) CASE "pllavfini": pllavfini! = VAL(val$): dfLOav! = pllavfini!'1p5 CASE "ffactor": ffactor! = VAL(val$)'1p5 CASE "srloavfac": srLOavfac! = VAL(val$)'1p5 CASE "srloavini": srLOavini! = VAL(val$)'1p5 END SELECT CASE "time" SELECT CASE name$ CASE "shortintegration": shortintegrationtime! = VAL(val$) CASE "longintegration": longintegrationfactor = VAL(val$) CASE "displaysamples": displaysamples = VAL(val$) CASE "utcoffset": utoffsh = VAL(val$) CASE "clockoffset": clockoffset! = VAL(val$) CASE "timercalavfac": timercalavfac! = val(val$) CASE "sleeptime": sleeptime! = VAL(val$) CASE "asf": asf! = VAL(val$) CASE "gri5990step": gri5990step# = VAL(val$) '1.043 milliseconds CASE "gri6731step": gri6731step# = VAL(val$) '0.28 milliseconds END SELECT CASE "amplitude" SELECT CASE name$ CASE "colour": colourflag = VAL(val$) CASE "gamma": gamma! = VAL(val$) CASE "dbrange": dbrange! = VAL(val$) CASE "fullscale" fullscale! = VAL(val$) IF RIGHT$(val$, 1) = "+" THEN linlogflag = 1 CASE "normalizedistance": normalizedistance! = VAL(val$) CASE "normalizepower": normalizepower! = VAL(val$) CASE "blankingoffsetcorr": blankingoffsetcorr! = VAL(val$) CASE "hipass": hipass! = VAL(val$) CASE "edgelevel": edgelevel! = VAL(val$) IF edgelevel! >= 1 THEN absthreshold! = edgelevel! CASE "phasereference": phasereference! = VAL(val$) CASE "referencestation": nstref = VAL(val$) CASE "firstblankedsample": firstblankedsample% = VAL(val$) CASE "lastblankedsample": lastblankedsample% = VAL(val$) CASE "notch100khz": notch100kHzflag% = VAL(val$) 'lcase khz case "xrcdelay7950": xrcdelay7950! = val(val$) '1.06 ms case "nb":nbflag%=val(val$) '0 NB off, 1 NB on, 2 only after Loran blanking case "nbinit":nbinit!=val(val$)^2 ' nb absolute threshold (lsb) case "nbfac":nbfac!=val(val$) 'nb power factor above average case "nbupdate": nbupdate!=val(val$) 'nb adaption (0..1) case "nbspread": nbspread%=val(val$) 'nb threshold to neighbour bins (0..2 left and right) case "nbcorr": nbcorr!=val(val$) 'nb intensity correction (0..1) END SELECT CASE "firfilter" SELECT CASE name$ CASE "hipass": hipass! = VAL(val$) CASE "firbegin": nfirbegin = VAL(val$) CASE "tap" afirr!(nfir) = VAL(val$) comma = INSTR(val$, ",") IF comma > 0 THEN afiri!(nfir) = VAL(RIGHT$(val$, LEN(val$) - comma)) 'PRINT afirr!(nfir), afiri!(nfir) nfir = nfir + 1 END SELECT CASE "taisync" SELECT CASE name$ CASE "searchrange": searchrange! = VAL(val$) CASE "integrationtime": avendtime# = VAL(val$) CASE "resynclimit": resynclimit! = VAL(val$) CASE "clockerrorlimit": clockerrorlimit! = VAL(val$) CASE "syncverbose": syncverbose% = 1 CASE "stations": naqm = VAL(val$) CASE "remove8000frompps":remove8000frompps%=1 CASE ELSE mxyz$ = UCASE$(RIGHT$(name$, 1)) griaq = VAL(name$) 'IF griaq > 4000 AND griag < 10000 AND INSTR("MVWXYZS", mxyz$) > 0 THEN IF griaq > 0 THEN griaq1A(nstaq) = griaq: mxyzaq1$(nstaq) = mxyz$ weightaq1!(nstaq) = VAL(val$) 'PRINT griaq1A(nstaq); mxyzaq1$(nstaq); weightaq1!(nstaq) nstaq = nstaq + 1 END IF END SELECT CASE "chaindelay" gried = VAL(name$) IF gried > 4000 AND gried < 10000 THEN gried1(ned) = gried ed1!(ned) = VAL(val$) 'given in milliseconds ned = ned + 1 END IF CASE "stations" SELECT CASE name$ CASE "hidden": hidden = VAL(val$) CASE "blankedstns": blankedstns = VAL(val$) CASE "blankingdistance": blankingdistance! = VAL(val$) CASE "singlepulses": singlepulses = VAL(val$) CASE "printstationlist": printstationlistflag% = 1 CASE ELSE 'mxyz$ = UCASE$(RIGHT$(name$, 1)) mf$=ltrim$(name$, ANY "-0123456789. ") 'remove number mxyz$=UCASE$(left$(mf$,1)) 'first character mf$=right$(mf$,len(mf$)-1) 'flags gri = VAL(name$): absgri = ABS(gri) IF absgri > 0 THEN gri1A(nst) = gri: mxyz1$(nst) = mxyz$: stnflags$(nst)=mf$ nst = nst + 1 END IF END SELECT END SELECT RETURN getblockofdata: IF INKEY$ = "q" THEN CLOSE 1: END IF actualsamples > 0 THEN ' > bufoffs THEN 'get average bufaverager& = 0&: bufaveragei& = 0& if nbflag% then bufaveragep!=0 FOR n% = 0 TO actualsamples - 1 bufaveragep! = bufaveragep! + apow!(n%) NEXT n% bufnbthre!=nbfac!*bufaveragep!/actualsamples FOR n% = 0 TO actualsamples - 1 if apow!(n%) inblockbytes THEN inbytes = inblockbytes ELSE inbytes = bytesleft in$ = INPUT$(inbytes, 1) bytesleft = bytesleft - inbytes actualsamples = LEN(in$) \ channels% \ bytespersample% 'flatter if-structures fixed out of heap space: IF channels% = 2 AND bytespersample% = 2 THEN FOR n = 0 TO actualsamples - 1 ar%(n + bufoffs) = CVI(MID$(in$, 1 + 4 * n, 2)) ai%(n + bufoffs) = CVI(MID$(in$, 3 + 4 * n, 2)) NEXT n END IF IF channels% = 2 AND bytespersample% = 1 THEN FOR n = 0 TO actualsamples - 1 ar%(n + bufoffs) = ASC(MID$(in$, 1 + 2 * n, 1)) - 128 ai%(n + bufoffs) = ASC(MID$(in$, 2 + 2 * n, 1)) - 128 NEXT n END IF IF channels% = 1 AND bytespersample% = 2 THEN FOR n = 0 TO actualsamples - 1 ar%(n + bufoffs) = CVI(MID$(in$, 1 + 2 * n, 2)) ai%(n + bufoffs) = 0 NEXT n END IF IF channels% = 1 AND bytespersample% = 1 THEN FOR n = 0 TO actualsamples - 1 ar%(n + bufoffs) = ASC(MID$(in$, 1 + n, 1)) - 128 ai%(n + bufoffs) = 0 NEXT n END IF IF filescale! <> 1! AND lofreq# = 0 THEN FOR n = 0 TO actualsamples - 1 ar%(n + bufoffs) = filescale! * ar%(n + bufoffs) IF channels% = 2 THEN ai%(n + bufoffs) = filescale! * ai%(n + bufoffs) NEXT n END IF IF lofreq# <> 0 OR lodrift! <> 0 THEN phcntm% = 100 lophasepersample! = 2! * pi! * lofreq# / srnominal! wwr! = COS(lophasepersample!): wwi! = -SIN(lophasepersample!) FOR n% = bufoffs TO bufoffs + actualsamples - 1 phcnt% = phcnt% + 1 IF phcnt% >= phcntm% THEN phcnt% = 0 phase! = 2! * pi! * locycles! wr! = filescale! * COS(phase!): wi! = filescale! * -SIN(phase!) 'conj.compl. locycles! = locycles! + lofreq# / srnominal! * phcntm% locycles! = locycles! - INT(locycles!) IF lodrift! <> 0 THEN lofreq# = lofreq# + lodrift! * phcntm% / srnominal! lophasepersample! = 2! * pi! * lofreq# / srnominal! wwr! = COS(lophasepersample!): wwi! = -SIN(lophasepersample!) END IF END IF IF channels% = 2 THEN ar! = ar%(n%): ai! = ai%(n%) ar%(n%) = ar! * wr! - ai! * wi! ai%(n%) = ar! * wi! + ai! * wr! ELSE ar! = ar%(n%) ar%(n%) = ar! * wr! ai%(n%) = ar! * wi! END IF wt! = wr! * wwr! - wi! * wwi! wi! = wr! * wwi! + wi! * wwr! wr! = wt! NEXT n% END IF if nbflag% then FOR n = bufoffs TO bufoffs + actualsamples - 1 ar!=ar%(n): ai!=ai%(n): apow!(n) = ar!*ar!+ai!*ai! next n end if bufstart# = bufend# - bufoffs / sr# bufend# = bufstart# + (bufoffs + actualsamples) / sr# srtime! = bufend# + starttime# - INT(starttime# / 86400) * 86400 'pctime! = killsec! - bytesleft / sr# / 4! - 3600 * utoffsh + 22 + leapsecs2k + clockoffset! '1q6: pctime! = killsec! - bytesleft / sr# / (channels% * bytespersample%) - 3600 * utoffsh + 22 + leapsecs2k + clockoffset! IF pctime! < 0 THEN pctime! = pctime! + 86400 pcsrdays = (pctime! - srtime!) / 86400 pcsrtime! = pctime! - srtime! - pcsrdays * 86400 'LOCATE , 1: PRINT USING "###.### "; srtime! - pctime!; IF bytesleft = 0 AND wavflag = 0 THEN CLOSE 1 'kill file KILL file1$ gosub gettimer if timercalcnt = 0 then gosub calibratetimer oldkillsec! = killsec! killsec!=todsec# -.05 GOSUB openfile END IF RETURN gettimer: timersec#=TIMER timepc$ = TIME$ todsec1&= val(left$(timepc$,2))*3600# + val(mid$(timepc$,4,2))*60# + val(right$(timepc$,2)) timercal# = timersec# - todsec1& - .5 if timercal#-timercalav# > 43200 then timercalav# = timercalav#+86400 if abs(timercal#-timercalav#) > .6 then if timercalcnt > 0 then timercalcnt = timercalcnt - 1 else timercalav# = timercalav# + (timercal#-timercalav#)*timercalavfac! if timercalcnt < 3 then timercalcnt = timercalcnt + 1 end if todsec# = timersec# - timercalav# '?timepc$;timersec#; timercal#;timercalav#;todsec# return calibratetimer: gosub gettimer oldtodsec1& = todsec1& while todsec1& = oldtodsec1& gosub gettimer sleep 50,1 'milliseconds wend timercalav# = timercal# +.5-.025 timercalcnt = 3 return openfile: IF wavflag = 0 THEN ON ERROR GOTO rxerr 'out of heap space appeared here OPEN file1$ FOR INPUT AS 1 ON ERROR GOTO 0 CLOSE 1 OPEN file1$ FOR BINARY AS 1 bytesleft = LOF(1) datduration! = bytesleft / sr# / 4 errordur! = datduration! - (killsec! - oldkillsec!) IF ABS(errordur!) < 2 THEN errorsum! = errorsum! + errordur! ' PRINT : PRINT killsec!; killsec! - oldkillsec!; datduration!; errordur!; errorsum! IF wavflag THEN hdrlen = 46 hdr$ = INPUT$(hdrlen, 1) bytesleft = bytesleft - hdrlen END IF RETURN rxerr: IF INKEY$ = "q" THEN END 'IF sleeptime! > 0 THEN SHELL "sleep " + STR$(sleeptime!) IF sleeptime! > 0 THEN SLEEP sleeptime! 'seconds RESUME printtime: secondsnow# = starttime# + bufend# - 22# - leapsecs2k - clockoffset! '? t100& = (secondsnow# - 86400# * INT(secondsnow# / 86400#)) * 100 hours% = t100& \ 360000 minutes% = (t100& - 360000 * hours%) \ 6000 seconds100% = (t100& - 360000 * hours% - 6000& * minutes%) hms100& = 1000000 * hours% + 10000& * minutes% + seconds100% hms$ = STR$(100000000 + hms100&) hms2$ = MID$(hms$, 3, 2) + ":" + MID$(hms$, 5, 2) + ":" + MID$(hms$, 7, 2)' + "." + MID$(hms$, 9, 2) PRINT hms2$; IF wavflag = 0 THEN PRINT USING " PC###.##s"; -pcsrtime!; 'PC clock offset hoursnow = hours% RETURN dataout: PRINT INT(integrationstart! + .5); "seconds" n = 0 FOR nst = 0 TO nstm - 1 PRINT "GRI"; griA(nst); "Stn"; mxyzA$(nst); avcntA(nst); "av" FOR nav = 0 TO displaysamples - 1 avr! = avr(n) / avcntA(nst) avi! = avi(n) / avcntA(nst) amag = SQR(avr! * avr! + avi! * avi!) PRINT amag; n = n + 1 NEXT nav PRINT NEXT nst RETURN writebmpheader: IF colourflag THEN depthbf = 3: colmaplen = 0 ELSE depthbf = 1: colmaplen = 256 * 4 END IF widt = nstdm * displaysamples widt = 4 * ((widt + 3) \ 4) 'must be multiple of 4 height = 0 'will be changed (offset 2 and 22) totallen = 54 + colmaplen + widt * height * depthbf hdr$ = "BM" + MKL$(totallen) + MKL$(0) + MKL$(54 + colmaplen) + MKL$(40) '18 hdr$ = hdr$ + MKL$(widt) + MKL$(height) + MKI$(1) + MKI$(depthbf * 8) '12 hdr$ = hdr$ + STRING$(24, 0) FOR ncol = 0 TO (colmaplen \ 4) - 1 hdr$ = hdr$ + (STRING$(3, ncol) + CHR$(0)) NEXT ncol OPEN file2$ FOR OUTPUT AS 2 PRINT #2, hdr$; CLOSE 2 bmpline$ = SPACE$(widt * depthbf) RETURN writebmpline: sqr3! = SQR(3!) n = 0 hipassflag% = ((hipass! <> 0) OR (hipassshift! <> 0)) hipassshf! = .5 * hipassshift! / longintegrationfactor xrcoffs1 = xrcoffs FOR nst = 0 TO nstdm - 1 scalefactor! = 1! / fullscale! IF normalizedistance! <> 0! AND dist!(nst) > 0 THEN scalefactor! = scalefactor! * dist!(nst) / normalizedistance! IF normalizepower! <> 0! AND pwrA!(nst) > 0 THEN scalefactor! = scalefactor! * SQR(normalizepower! / pwrA!(nst)) hipass1! = -hipass! - hipassshf! * serrorl!(nst) hipass2! = -hipass! + hipassshf! * serrorl!(nst) '======= nxrcm=0 IF (griA(nst)=6000 AND xrccnt6000=3) THEN nxrcm=2 'IF (griA(nst)=7950 AND xrccnt7950=15) THEN nxrcm=4 IF (griA(nst)=7950 AND xrccnt7950=9) THEN nxrcm=2 IF xrcoffs1 < 0 then FOR nxrc= 0 to nxrcm-1 nxrcd = (nstm + xrcoffs1)* displaysamples FOR nav = 0 TO displaysamples - 1 avrA!(n + nav) = avrA!(n + nav) - avrA!(nxrcd + nav) / 2 aviA!(n + nav) = aviA!(n + nav) - aviA!(nxrcd + nav) / 2 NEXT nav xrcoffs1 = xrcoffs1 + 1 next nxrc end if '... '6000m '... '7950m '7950w '7950x '... 'hidden=14 '6000m '6000m '7950m '7950m '7950w '7950w '7950x '7950x '[end] '======== FOR nav = 0 TO displaysamples - 1 avr! = avrA!(n): avi! = aviA!(n) IF hipassflag% THEN IF nav > 0 AND nav < displaysamples - 1 THEN avr! = avr! + hipass1! * avrA!(n - 1) + hipass2! * avrA!(n + 1) avi! = avi! + hipass1! * aviA!(n - 1) + hipass2! * aviA!(n + 1) END IF END IF IF nfirm > 0 THEN ndfir1 = -nfirbegin IF ndfir1 + nav >= displaysamples THEN ndfir1 = displaysamples - 1 - nav ndfir2 = -(nfirbegin + nfirm - 1) IF ndfir2 + nav < 0 THEN ndfir2 = -nav FOR ndfir = ndfir1 TO ndfir2 STEP -1 nfir = -ndfir - nfirbegin avr! = avr! + avrA!(n + ndfir) * afirr!(nfir) - aviA!(n + ndfir) * afiri!(nfir) avi! = avi! + avrA!(n + ndfir) * afiri!(nfir) + aviA!(n + ndfir) * afirr!(nfir) NEXT ndfir END IF avr! = scalefactor! * avr! / longintegrationfactor avi! = scalefactor! * avi! / longintegrationfactor amag! = SQR(avr! * avr! + avi! * avi!) lum! = amag! IF dbrange! = 0! THEN IF gamma! <> 1! THEN lum! = lum! ^ gamma! IF linlogflag AND lum! > .75 THEN lum! = .75 + 8.6859 * LOG(lum! / .75) / 256! ELSE IF amag! > 0 THEN lum! = 1! + 8.6859 * LOG(lum!) / ABS(dbrange!) IF lum! < 0! THEN lum! = 0! END IF IF colourflag = 0 THEN IF lum! > 1! THEN lum! = 1! MID$(bmpline$, 1 + n, 1) = CHR$(INT(255.999 * lum!)) ELSE IF amag! > 0! THEN colr! = avr! / amag!: coli! = avi! / amag! colmag! = 1! IF colourflag = 1 THEN lum! = lum! / (lum! + 1!) ELSE lum! = .5 * lum!: IF lum! > .5 THEN lum! = .5 END IF col2! = 2! * colr! IF col2! > colmag! THEN col2! = colmag! IF col2! < -colmag! THEN col2! = -colmag! IF lum! > .5 THEN col2! = (1! / lum! - 1!) * col2! red = INT(255.999 * lum! * (1! + col2!)) col2! = sqr3! * coli! - colr! IF col2! > colmag! THEN col2! = colmag! IF col2! < -colmag! THEN col2! = -colmag! IF lum! > .5 THEN col2! = (1! / lum! - 1!) * col2! blue = INT(255.999 * lum! * (1! + col2!)) col2! = -sqr3! * coli! - colr! IF col2! > colmag! THEN col2! = colmag! IF col2! < -colmag! THEN col2! = -colmag! IF lum! > .5 THEN col2! = (1! / lum! - 1!) * col2! green = INT(255.999 * lum! * (1! + col2!)) bgr$ = CHR$(blue) + CHR$(green) + CHR$(red) MID$(bmpline$, 1 + 3 * n, 3) = bgr$ END IF n = n + 1 NEXT nav NEXT nst newhourflag%=0 IF hoursnow <> hoursnowold THEN newhourflag%=1 hoursnowold = hoursnow tickb = depthbf * (1 - (hoursnow = 0)) '2 pix at midnight FOR nst = 0 TO nstdm - 1 MID$(bmpline$, 1 + nst * displaysamples * depthbf, tickb) = STRING$(tickb, 255) NEXT nst if hoursnow=0 and gri5990step#<>0 then 'gri5990 step at midnight for nst = 0 TO nstm - 1 if griA(nst)=5990 then nextpulsestartA#(nst) = nextpulsestartA#(nst) + gri5990step# * 0.001# firstpulsestartA#(nst) = firstpulsestartA#(nst) + gri5990step# * 0.001# '**** 1q17 txdelayA!(nst)=txdelayA!(nst) + gri5990step# 'txdelayaqA!(nst)=txdelayaqA!(nst) + gri5990step# 'removed in 1q18 end if next nst end if if hoursnow=0 and gri6731step#<>0 then 'gri6731 step at midnight for nst = 0 TO nstm - 1 if griA(nst)=6731 then nextpulsestartA#(nst) = nextpulsestartA#(nst) + gri6731step# * 0.001# firstpulsestartA#(nst) = firstpulsestartA#(nst) + gri6731step# * 0.001# txdelayA!(nst)=txdelayA!(nst) + gri6731step# end if next nst end if END IF totallen = totallen + depthbf * widt height = height + 1 bmpretry: OPEN file2$ FOR BINARY AS 2 if bmperrorhandling% then ON ERROR GOTO bmperr PUT #2, 3, totallen '<=1q3 used to crash here PUT #2, 23, height PUT #2, totallen - depthbf * widt + 1, bmpline$ ON ERROR GOTO 0 '1q4 CLOSE 2 if bmpcmd$<>"" then if (bmpcmdcount>0 or newhourflag%>0) and bmpcnt>0 then bmpcnt = bmpcnt - 1 if bmpcnt=0 then bmpcnt=abs(bmpcmdcount) shell bmpcmd$ end if end if RETURN bmperr: 'this keeps LV running but may occasionally cause loss of bmp contents! close 2 IF INKEY$ = "q" THEN END if bmperrorhandling% > 1 then print err sleep 200,1 '200ms goto bmpretry measureslope: threshold! = edgelevel! '.45 GOSUB getedgedata t1cyc! = ((t1! - centralsample + 1) / sr# + dtexp2#(nstref)) * 100000 phasecyc! = phase! / 2! / pi! IF avcnt! > 0 THEN serrorcyc! = serrorA!(nstref) / avcnt! / sr# * 100000! '+-10us for 7499 22042Hz sr ELSE serrorcyc! = 0 END IF IF srcorrfac! THEN deltasr1# = srcorrfac! * (t1! - centralsample + 1) / shortintegrationtime! deltabuftime# = -srcorrpfac! * (t1! - centralsample + 1)/sr# bufstart# = bufstart# + deltabuftime# bufend# = bufend# + deltabuftime# sr# = sr# + deltasr1# END IF dflo! = phasecyc! - phasecycold! dflo! = (dflo! - INT(dflo! + .5)) / shortintegrationtime! / 100000! phasecnt = phasecnt - INT(phasecyc! - phasecycold! + .5) phasecycold! = phasecyc!: LOus! = LOus! + dflo! * shortintegrationtime! * 1000000! 'correct LO phasetime for sr deviation srLOus! = ((sr# - sr0#) * srlofactor! + srlooffset! * 1000000!) * shortintegrationtime! LOus! = LOus! - srLOus! * 1 IF amagmax! > resynclimit! THEN validcnt = validcnt + 1 ELSE validcnt = 0 IF validcnt >= 4 THEN validcnt = 4: lastLOus! = LOus! ELSE LOus! = 0: srLOav! = srLOavini! END IF IF validcnt = 4 AND pllflag% THEN dfLOav! = dflo! * pllavf! + dfLOav! * (1 - pllavf!) 'corrected 1p5 pllf# = pllcenter# - pllfactort! * LOus! - pllfactorf! * (dflo! - dfLOav!) IF srpllfactor! THEN pllf# = pllf# * (1 + srpllfactor! * (sr0# / sr# - 1)) srLOav! = srLOavfac! * srLOus! / shortintegrationtime! / 1000000 + (1 - srLOavfac!) * srLOav! '1p5 IF ffactor! THEN pllf# = pllf# + ffactor! * (dfLOav! - srLOav!)'1p5 IF plldither! THEN pllf# = pllf# + plldither! * (2 * RND - 1) IF pllf# > pllhigh! THEN pllf# = pllhigh! IF pllf# < plllow! THEN pllf# = plllow! plltestsig$ = " " + LEFT$(STR$(pllf#), 10) SHELL pllcmd$ + plltestsig$ END IF PRINT USING " max##### edge##.## "; amagmax!; t1!; PRINT USING "SR #####.####Hz "; sr#; PRINT USING "LO +#.#^^^^######.##us"; dflo!; lastLOus!; LOCATE CSRLIN-1,1 : PRINT '1q5 IF phasecorrflag THEN IF amagslope! > 0 THEN phasefacr! = (asloper! * phaserefr! - aslopei! * phaserefi!) / amagslope! phasefaci! = (aslopei! * phaserefr! + asloper! * phaserefi!) / amagslope! END IF END IF IF notch100kHzflag% THEN IF bufaveragen! > 0 THEN offset100kr! = bufaverager! / bufaveragen! offset100ki! = bufaveragei! / bufaveragen! END IF bufaverager! = 0: bufaveragei! = 0: bufaveragen! = 0 END IF offsetrb! = offsetr!: offsetib! = offseti! FOR nst = 0 TO nstm - 1 avcnt! = avcntA(nst) IF avcnt! > 0 THEN avcntbl! = avcnt!* blankingfactorA!(nst) serrorl!(nst) = serrorl!(nst) + serrorA!(nst) / avcnt! IF notch100kHzflag% THEN factor! = avcntsignA(nst) / avcntbl! offsetrb4! = offset100kr! * factor! offsetib4! = offset100ki! * factor! ELSE offsetrb4! = offsetrb! * 4! / npulsemA(nst) offsetib4! = offsetib! * 4! / npulsemA(nst) END IF '?nst;avcnt!; FOR n = nst * displaysamples TO (nst + 1) * displaysamples - 1 avr1! = avr(n) / avcntbl! - offsetrb4! avi1! = avi(n) / avcntbl! - offsetib4! if nbflag%>0 and nbcorr!>0 then factor!=(1! - nbcorr!*nbcnt!(n)/avcnt!) nbcnt!(n)=0! '?factor!; if factor!>0 then avr1!=avr1!/factor! avi1!=avi1!/factor! end if end if avrA!(n) = avrA!(n) + avr1! * phasefacr! + avi1! * phasefaci! aviA!(n) = aviA!(n) + avi1! * phasefacr! - avr1! * phasefaci! NEXT n '? IF blankingflag%(nst) THEN '4 unbalanced pulses per 2xGRI avblr! = 0!: avbli! = 0! FOR n = nst * displaysamples + firstblankedsample% TO nst * displaysamples + lastblankedsample% avblr! = avblr! + avr(n): avbli! = avbli! + avi(n) NEXT n blankedsampleoffsfac! = blankingoffsetcorr! * 4! / sr# / griA(nst) * 50000! / avcnt! offsetrb! = offsetrb! - avblr! * blankedsampleoffsfac! offsetib! = offsetib! - avbli! * blankedsampleoffsfac! END IF END IF avcntA(nst) = 0 avcntsignA(nst) = 0 serrorA!(nst) = 0! NEXT nst FOR n = 0 TO nstm * displaysamples - 1 avr(n) = 0: avi(n) = 0 NEXT n IF amagmax! < resynclimit! OR (wavflag = 0 AND clockerrorlimit! > 0 AND ABS(pcsrtime!) > clockerrorlimit!) THEN resynccnt = resynccnt + 1 ELSE IF resynccnt > 0 THEN resynccnt = resynccnt - 1 END IF RETURN getedgedata: 'input: nstref, threshold! 'output: amagmax!,t1!,amagslope!,asloper!,aslopei!,phase!, avcnt! nav0 = nstref * displaysamples avcnt! = avcntA(nstref) IF avcnt! = 0! THEN RETURN nav1 = centralsample - 4: IF nav1 < 0 THEN nav1 = 0 nav2 = centralsample + 4: IF nav2 >= displaysamples THEN nav2 = displaysamples - 1 amagmax! = 0: nmax = centralsample FOR nav = nav1 TO nav2 avr! = avr(nav0 + nav) / avcnt!: avi! = avi(nav0 + nav) / avcnt! amag! = SQR(avr! * avr! + avi! * avi!) amagA!(nav) = amag! IF amag! > amagmax! THEN amagmax! = amag!: nmax = nav NEXT nav slopelim! = amagmax! * threshold! IF absthreshold! > 0 THEN slopelim! = absthreshold! 'n1 = nmax - 1: amagA!(0) = 0! 'backward search from peak 'WHILE amagA!(n1) > slopelim!: n1 = n1 - 1: WEND n1 = 1: amagA!(0) = 0! 'forward search WHILE amagA!(n1) < slopelim! AND n1 < nmax: n1 = n1 + 1: WEND n1 = n1 - 1 IF amagA!(n1 + 1) > amagA!(n1) THEN f1! = (slopelim! - amagA!(n1)) / (amagA!(n1 + 1) - amagA!(n1)) IF f1! < 0 OR f1! > 1 THEN f1! = .5 f0! = 1 - f1! asloper! = (f0! * avr(nav0 + n1) + f1! * avr(nav0 + n1 + 1)) / avcnt! aslopei! = (f0! * avi(nav0 + n1) + f1! * avi(nav0 + n1 + 1)) / avcnt! t1! = n1 + f1! - serrorA!(nstref) / avcnt! amagslope! = SQR(asloper! * asloper! + aslopei! * aslopei!) IF amagslope! > 0 THEN IF ABS(asloper!) > ABS(aslopei!) THEN phase! = ATN(aslopei! / asloper!) ELSE phase! = .5 * pi! - ATN(asloper! / aslopei!) END IF IF (asloper! + aslopei!) < 0 THEN phase! = phase! + pi! IF phase! < 0 THEN phase! = phase! + 2! * pi! END IF RETURN printstatistics: dtmean# = sumdt1# / sumdt0# dtmean! = dtmean# dtstdev! = SQR(ABS(sumdt2# / sumdt0# - dtmean# * dtmean#)) PRINT "envelope diff mean"; dtmean!; "st dev"; dtstdev!; "rms" RETURN synchronize: IF avendtime# = 0 THEN RETURN aqbufend = 0 FOR nst = 0 TO nstaqm - 1 aqbufstartA(nst) = aqbufend aqbufptrA(nst) = aqbufend aqbufend = aqbufend + griaqA(nst) / 50000# * sr# aqbufendA(nst) = aqbufend bufptr(nst) = bufoffs startaq#(nst) = bufoffs / sr# NEXT nst DIM avaqr(aqbufend - 1), avaqi(aqbufend - 1) PRINT PRINT "sync: averaging for"; INT(avendtime# + .5); "sec "; LOCATE , 1 totalsamples = 0 FOR nst = 0 TO nstaqm - 1 avcntaq(nst) = 0 NEXT nst WHILE bufend# < avendtime# GOSUB getblockofdata FOR nst = 0 TO nstaqm - 1 aqbufend = aqbufendA(nst) naq = aqbufptrA(nst) IF bandcenter# THEN firstsample = totalsamples + bufptr(nst) - bufoffs phase# = firstsample / sr# * looffsetaq#(nst) phase! = 2 * pi# * (phase# - INT(phase#)) omega! = 2 * pi# * looffsetaq#(nst) / sr# END IF FOR nbuf = bufptr(nst) TO bufoffs + actualsamples - 1 IF bandcenter# THEN cosp! = COS(phase!): sinp! = -SIN(phase!) ar! = ar%(nbuf): ai! = ai%(nbuf) ar = ar! * cosp! - ai! * sinp! ai = ar! * sinp! + ai! * cosp! avaqr(naq) = avaqr(naq) + ar avaqi(naq) = avaqi(naq) + ai phase! = phase! + omega! IF phase! > pi! THEN phase! = phase! - (2 * pi!) 'IF nst = 0 THEN PRINT phase!; ELSE avaqr(naq) = avaqr(naq) + ar%(nbuf) avaqi(naq) = avaqi(naq) + ai%(nbuf) END IF naq = naq + 1 IF naq >= aqbufend THEN startaq#(nst) = startaq#(nst) + griaqA(nst) / 50000# buftime# = bufstart# + (bufoffs + nbuf) / sr# delta = (startaq#(nst) - buftime#) * sr# '+-1 samples nbuf = nbuf + delta naq = aqbufstartA(nst) avcntaq(nst) = avcntaq(nst) + 1 END IF NEXT nbuf bufptr(nst) = nbuf - actualsamples 'bufoffs-+1 if last delta was -+1 aqbufptrA(nst) = naq NEXT nst totalsamples = totalsamples + actualsamples WEND 'correlate to pulsegroups aqbufend = aqbufendA(nstaqm - 1) FOR nst = 0 TO nstaqm - 1 FOR ms = 0 TO 1 '0 master, 1 secondary FOR np = 0 TO 7 pulseoffscor(nst, ms, np) = 100 * np / 100000# * sr# NEXT np FOR np = 0 TO 7 pulseoffscor(nst, ms, 9 - ms + np) = (100 * np + griaqA(nst)) / 100000# * sr# NEXT np IF ms = 0 THEN pulseoffscor(nst, 0, 8) = 900 / 100000# * sr# pulseoffscor(nst, 0, 17) = (900 + griaqA(nst)) / 100000# * sr# END IF NEXT ms FOR np = 0 TO 17 pulsesigncor%(nst, 0, np) = 2 * VAL(MID$(searchpattern$, 1 + np, 1)) - 1 NEXT np FOR np = 0 TO 15 pulsesigncor%(nst, 1, np) = 2 * VAL(MID$(searchpattern$, 19 + np, 1)) - 1 NEXT np NEXT nst maxgriaq = 0 FOR nst = 0 TO nstaqm - 1 IF griaqA(nst) > maxgriaq AND UCASE$(mxyzaq$(nst)) <> "A" THEN maxgriaq = griaqA(nst): maxnst = nst NEXT nst corbufend = aqbufendA(maxnst) - aqbufstartA(maxnst) DIM cormpA!(corbufend - 1), corsp!(corbufend - 1) FOR nst = 0 TO nstaqm - 1 aqbufstart = aqbufstartA(nst) aqbufend = aqbufendA(nst) corend = aqbufend - aqbufstart IF mxyzaq$(nst) = "A" THEN 'single pulse maxpwr! = 0 naq = aqbufstart FOR nc = 0 TO corend - 1 cmr! = avaqr(aqbufstart + nc): cmi! = avaqi(aqbufstart + nc) if griaqA(nst)=50000 and remove8000frompps% then nc1=nc-sr#*.04: if nc1<0 then nc1=nc1+corend cmr! = cmr!-.5*avaqr(aqbufstart + nc1): cmi! = cmi!-.5*avaqi(aqbufstart + nc1) nc1=nc+sr#*.04: if nc1>=corend then nc1=nc1-corend cmr! = cmr!-.5*avaqr(aqbufstart + nc1): cmi! = cmi!-.5*avaqi(aqbufstart + nc1) end if cormp! = cmr! * cmr! + cmi! * cmi! IF cormp! > maxpwr! THEN maxpwr! = cormp!: maxnc = nc NEXT nc ampl! = SQR(maxpwr!) / avcntaq(nst): delay! = maxnc / sr# amplA!(nst, 0) = ampl!: delayA!(nst, 0) = delay! IF syncverbose% THEN PRINT griaqA(nst); "max"; ampl!; "at"; delay! * 1000; "ms" END IF ELSE 'pulsegroups FOR nc = 0 TO corend - 1 cmr! = 0: cmi! = 0 FOR np = 0 TO 17 naq = aqbufstart + nc + pulseoffscor(nst, 0, np) IF naq >= aqbufend THEN naq = naq - corend 'aqbufstart IF pulsesigncor%(nst, 0, np) > 0 THEN cmr! = cmr! + avaqr(naq) cmi! = cmi! + avaqi(naq) ELSE cmr! = cmr! - avaqr(naq) cmi! = cmi! - avaqi(naq) END IF NEXT np cormpA!(nc) = (cmr! * cmr! + cmi! * cmi!) / 324! cmr! = 0: cmi! = 0 FOR np = 0 TO 15 naq = aqbufstart + nc + pulseoffscor(nst, 1, np) IF naq >= aqbufend THEN naq = naq - corend 'aqbufstart IF pulsesigncor%(nst, 1, np) > 0 THEN cmr! = cmr! + avaqr(naq) cmi! = cmi! + avaqi(naq) ELSE cmr! = cmr! - avaqr(naq) cmi! = cmi! - avaqi(naq) END IF NEXT np corsp!(nc) = (cmr! * cmr! + cmi! * cmi!) / 256! NEXT nc mcnt = 0: scnt = 0 'masters and slaves found 'delete old excess slaves FOR nmsss = 0 TO nmsssm - 1 amplA!(nst, nmsss) = 0 NEXT nmsss 'search peaks, strongest first FOR nmsss = 0 TO nmsssm - 1 maxpwr! = 0 FOR nc = 0 TO corend - 1 IF cormpA!(nc) > maxpwr! THEN maxpwr! = cormpA!(nc): maxnc = nc: maxms = 0 IF corsp!(nc) > maxpwr! THEN maxpwr! = corsp!(nc): maxnc = nc: maxms = 1 NEXT nc IF maxms = 0 THEN mssi = 0: mcnt = mcnt + 1 IF mcnt > 1 THEN mssi = 100 ELSE scnt = scnt + 1: mssi = scnt END IF ampl! = SQR(maxpwr!) / avcntaq(nst): delay! = maxnc / sr# IF mssi < nmsssm THEN amplA!(nst, mssi) = ampl!: delayA!(nst, mssi) = delay! END IF 'clear sidelobes +-9ms clearw = .009 * sr# nc = maxnc - clearw: IF nc < 0 THEN nc = nc + corend FOR ncc = 0 TO 2 * clearw cormpA!(nc) = 0: corsp!(nc) = 0 nc = nc + 1: IF nc >= corend THEN nc = 0 NEXT ncc nc = maxnc + (corend / 2) - clearw: IF nc >= corend THEN nc = nc - corend FOR ncc = 0 TO 2 * clearw cormpA!(nc) = 0: corsp!(nc) = 0 nc = nc + 1: IF nc >= corend THEN nc = 0 NEXT ncc IF syncverbose% THEN IF maxms THEN PRINT "S"; ELSE PRINT "M"; PRINT griaqA(nst); "max"; ampl!; "at"; delay! * 1000; "ms" END IF 'IF nst = 2 AND nmsss = 0 THEN GOSUB outfiles NEXT nmsss END IF NEXT nst 'calctai: assec! = assumedstart# - 86400 * INT(assumedstart# / 86400) pcsec! = pctime! - totalsamples / sr# IF wavflag = 0 THEN assumedstart# = assumedstart# + pcsec! - assec! + .05 maxerror! = searchrange! FOR nst = 0 TO nstaqm - 1 t2gri# = griaqA(nst) / 50000# deldraq!(nst) = txdelayaqA!(nst) / 1000! + distaq!(nst) * 1000! / c! expdelgris# = (deldraq!(nst) - assumedstart#) / t2gri# dtcalcaq!(nst) = (expdelgris# - INT(expdelgris#)) * t2gri# IF INSTR("AKLMN", UCASE$(mxyzaq$(nst))) > 0 THEN nmsss = 0 ELSE nmsss = 1 dtmeasaq!(nst) = delayA!(nst, nmsss) NEXT nst t2gri0# = griaqA(0) / 50000# startoffs# = -dtmeasaq!(0) + dtcalcaq!(0) if startoffs# > t2gri0# / 2 then startoffs# = startoffs# - t2gri0# '1q7 if startoffs# < -t2gri0# / 2 then startoffs# = startoffs# + t2gri0# '1q7 ndt1m = maxerror! / t2gri0# mindelta! = 10 ^ 6 deltalim! = 1 FOR ndt1 = 0 TO ndt1m FOR ndtsign = 0 TO 1 ndt = ndt1 * (2 * ndtsign - 1) startoffsdt# = startoffs# + ndt * t2gri0# sumqdelta! = 0 FOR nst = 1 TO nstaqm - 1 t2gri# = griaqA(nst) / 50000# deltagris# = (dtcalcaq!(nst) - dtmeasaq!(nst) - startoffsdt#) / t2gri# delta! = (deltagris# - INT(deltagris# + .5#)) * t2gri# 'PRINT griaqA(nst); ndt; "delta"; delta! * 1000! sumqdelta! = sumqdelta! + delta! * delta! * weightaq!(nst) * weightaq!(nst) NEXT nst IF sumqdelta! < mindelta! THEN mindelta! = sumqdelta!: ndtmin = ndt IF sumqdelta! < deltalim! THEN ndtsign = 1: ndt1 = ndt1m NEXT ndtsign NEXT ndt1 startoffsmin# = startoffs# IF mindelta! < deltalim! THEN startoffsmin# = startoffs# + ndtmin * t2gri0# '1q8 LOCATE , 1 IF mindelta! < deltalim! THEN PRINT "TAI match"; ELSE PRINT "TAI guess"; PRINT USING " at ####.## s offset, rel. error ##.## rms, "; startoffsmin#; SQR(mindelta!); starttime# = assumedstart# + startoffsmin# startday# = (starttime# - 22 - leapsecs2k - clockoffset!) / 86400# startday = INT(startday#) starthour# = (startday# - startday) * 24: starthour = INT(starthour#) startmin# = (starthour# - starthour) * 60: startmin = INT(startmin#) startsec! = (startmin# - startmin) * 60 s1000& = 1000 * startsec! hms1000& = 10000000 * starthour + 100000 * startmin + s1000& hms$ = STR$(1000000000 + hms1000&) hms2$ = MID$(hms$, 3, 2) + ":" + MID$(hms$, 5, 2) + ":" + MID$(hms$, 7, 2) + "." + MID$(hms$, 9, 3) PRINT " "; hms2$; PRINT " UT" ERASE avaqr, avaqi, cormpA!, corsp! RETURN display: PRINT griaqA(nst); avcntaq(nst); "av" FOR n = aqbufstartA(nst) TO aqbufendA(nst) - 1 'FOR n = 0 TO aqbufendA(nst) - aqbufstartA(nst) - 1 ar! = avaqr(n) / avcntaq(nst): ai! = avaqi(n) / avcntaq(nst) 'ar! = cormr(n) / avcntaq(nst) / 18: ai! = cormi(n) / avcntaq(nst) / 18 'ar! = corsr(n) / avcntaq(nst): ai! = corsi(n) / avcntaq(nst) aqmag = SQR(ar! * ar! + ai! * ai!) * 1 aqmag$ = "_" IF aqmag > 50 THEN aqmag$ = "-" IF aqmag > 200 THEN aqmag$ = "^" PRINT aqmag$; NEXT n PRINT RETURN outfiles: OPEN "avfile.txt" FOR OUTPUT AS 3 FOR n = aqbufstartA(nst) TO aqbufendA(nst) - 1 PRINT #3, avaqr(n); avaqi(n) NEXT n CLOSE 3 'OPEN "corfile.txt" FOR OUTPUT AS 3 'FOR n = 0 TO aqbufendA(nst) - aqbufstartA(nst) - 1 'PRINT #3, cormr(n); cormi(n); corsr(n); corsi(n) 'NEXT n 'CLOSE 3 OPEN "corpfile.txt" FOR OUTPUT AS 3 FOR n = 0 TO aqbufendA(nst) - aqbufstartA(nst) - 1 PRINT #3, cormpA!(n); corsp!(n) NEXT n CLOSE 3 RETURN evaluatedateandtime: 'date1$="mm-dd-yyyy", time1$="hh:mm:ss" datey = VAL(RIGHT$(date1$, 4)) datem = VAL(LEFT$(date1$, 2)) dated = VAL(MID$(date1$, 4, 2)) taidays = (datey - 1958) * 365 + 11 taidays = taidays + INT((datey - 2001) / 4) taidays = taidays - INT((datey - 2001) / 100) monthlen: DATA 31,28,31,30,31,30,31,31,30,31,30,31 RESTORE monthlen FOR month = 2 TO datem READ daysinmonth taidays = taidays + daysinmonth NEXT month leapsecs2k = 0 IF datey > 2005 THEN leapsecs2k = leapsecs2k + 1 IF datey > 2008 THEN leapsecs2k = leapsecs2k + 1 IF datey + datem / 12! >= 2012.5 THEN leapsecs2k = leapsecs2k + 1 inleapyear = (datey - 4 * INT(datey / 4) = 0) inleapyear = inleapyear AND (datey - 100 * INT(datey / 100) <> 0) IF inleapyear AND datem > 2 THEN taidays = taidays + 1 taidays = taidays + dated - 1 timeh = VAL(LEFT$(time1$, 2)) timem = VAL(MID$(time1$, 4, 2)) times = VAL(RIGHT$(time1$, 2)) taisecs! = (timeh - utoffsh) * 3600 + timem * 60 + times + 32 + leapsecs2k + clockoffset! assumedstart# = taidays * 3600# * 24# + taisecs! - 10 RETURN calculatedistance: 'input lat#,lon#,mylat#,mylon# (degrees) 'output dist#,alf#,bet# (km, degrees NESW) pi# = 4# * ATN(1) rearth# = 6370# gam# = (lon# - mylon#) * pi# / 180 a# = mylat# * pi# / 180#: b# = lat# * pi# / 180# cosc# = SIN(a#) * SIN(b#) + COS(a#) * COS(b#) * COS(gam#) sinc# = SQR(1# - cosc# * cosc#) cc# = ATN(sinc# / cosc#) IF cc# < 0 THEN cc# = cc# + pi# dist# = cc# * rearth# cosalf# = COS(a#) * COS(b#) * SIN(gam#) sinalf# = COS(a#) * SIN(b#) - SIN(a#) * COS(b#) * COS(gam#) IF ABS(cosalf#) > ABS(sinalf#) THEN alf# = ATN(sinalf# / cosalf#) IF cosalf# < 0 THEN alf# = alf# - pi# ELSE alf# = pi# / 2 - ATN(cosalf# / sinalf#) IF sinalf# < 0 THEN alf# = alf# + pi# END IF alf# = pi# / 2 - alf# alf# = alf# - 2 * pi# * INT(alf# / (2 * pi#)) alf# = alf# * 180# / pi# cosbet# = -cosalf# sinbet# = COS(b#) * SIN(a#) - SIN(b#) * COS(a#) * COS(gam#) IF ABS(cosbet#) > ABS(sinbet#) THEN bet# = ATN(sinbet# / cosbet#) IF cosbet# < 0 THEN bet# = bet# - pi# ELSE bet# = pi# / 2 - ATN(cosbet# / sinbet#) IF sinbet# < 0 THEN bet# = bet# + pi# END IF bet# = pi# / 2 - bet# bet# = bet# - 2 * pi# * INT(bet# / (2 * pi#)) bet# = bet# * 180# / pi# RETURN calculatevincenty: ' adapted from from chris veness ' http://www.movable-type.co.uk/scripts/latlong-vincenty.html Lat1# = mylat# Lon1# = mylon# Lat2# = lat# Lon2# = lon# a# = 6378137 f# = 1 / 298.257223563# b# = a# * (1# - f#) eps# = 1E-12 pi# = 4# * ATN(1#) L# = (Lon2# - Lon1#) * pi# / 180 U1# = ATN((1# - f#) * TAN(Lat1# * pi# / 180)) sinU1# = SIN(U1#) cosU1# = COS(U1#) U2# = ATN((1 - f#) * TAN(Lat2# * pi# / 180)) sinU2# = SIN(U2#) cosU2# = COS(U2#) lambda# = L# lambdaP# = 2 * pi# iterlimit% = 20 WHILE ABS(lambda# - lambdaP#) > eps# AND iterlimit% > 0 iterlimit% = iterlimit% - 1 sinlambda# = SIN(lambda#) coslambda# = COS(lambda#) x1# = cosU2# * sinlambda# x2# = cosU1# * sinU2# - sinU1# * cosU2# * coslambda# sinsigma# = SQR(x1# * x1# + x2# * x2#) IF sinsigma# = 0 THEN dist# = 0: alf# = 0: bet# = 180 RETURN END IF cossigma# = sinU1# * sinU2# + cosU1# * cosU2# * coslambda# IF ABS(cossigma#) > ABS(sinsigma#) THEN sigma# = ATN(sinsigma# / cossigma#) IF cossigma# < 0 AND sinsigma# < 0 THEN sigma# = sigma# - pi# IF cossigma# < 0 AND sinsigma# > 0 THEN sigma# = sigma# + pi# ELSE sigma# = pi# / 2 - ATN(cossigma# / sinsigma#) IF sinsigma# < 0 THEN sigma# = sigma# + pi# END IF sinalpha# = cosU1# * cosU2# * sinlambda# / sinsigma# cosSqalpha# = 1# - sinalpha# * sinalpha# IF cosSqalpha# <> 0# THEN cos2sigmaM# = cossigma# - 2 * sinU1# * sinU2# / cosSqalpha# ELSE cos2sigmaM# = 0 END IF c# = f# / 16 * cosSqalpha# * (4 + f# * (4 - 3 * cosSqalpha#)) lambdaP# = lambda# x# = cos2sigmaM# + c# * cossigma# * (-1# + 2 * cos2sigmaM# * cos2sigmaM#) lambda# = L# + (1# - c#) * f# * sinalpha# * (sigma# + c# * sinsigma# * x#) WEND uSq# = cosSqalpha# * (a# * a# / (b# * b#) - 1#) AA# = 1# + uSq# / 16384 * (4096 + uSq# * (-768 + uSq# * (320 - 175 * uSq#))) BB# = uSq# / 1024 * (256 + uSq# * (-128 + uSq# * (74 - 47 * uSq#))) x# = cos2sigmaM# * (4 * sinsigma# ^ 2 - 3) * (4 * cos2sigmaM# ^ 2 - 3) x# = BB# / 4 * (cossigma# * (2 * cos2sigmaM# ^ 2 - 1) - BB# / 6 * x#) deltasigma# = BB# * sinsigma# * (cos2sigmaM# + x#) s# = b# * AA# * (sigma# - deltasigma#) dist# = s# / 1000 sinalf# = cosU2# * sinlambda# cosalf# = cosU1# * sinU2# - sinU1# * cosU2# * coslambda# IF ABS(cosalf#) > ABS(sinalf#) THEN alf# = ATN(sinalf# / cosalf#) IF cosalf# < 0 THEN alf# = alf# - pi# ELSE alf# = pi# / 2 - ATN(cosalf# / sinalf#) IF sinalf# < 0 THEN alf# = alf# + pi# END IF alf# = alf# - 2 * pi# * INT(alf# / (2 * pi#) + m5) alf# = alf# * 180# / pi# sinbet# = cosU1# * sinlambda# cosbet# = -sinU1# * cosU2# + cosU1# * sinU2# * coslambda# IF ABS(cosbet#) > ABS(sinbet#) THEN bet# = ATN(sinbet# / cosbet#) IF cosbet# < 0 THEN bet# = bet# - pi# ELSE bet# = pi# / 2 - ATN(cosbet# / sinbet#) IF sinbet# < 0 THEN bet# = bet# + pi# END IF bet# = pi# + bet# bet# = bet# - 2 * pi# * INT(bet# / (2 * pi#) + .5) bet# = bet# * 180# / pi# RETURN