!************************************************************* ! XLDAM : Extract Line data from daq2asc output * ! with real-time Pos. data interporation * ! and calculation of IGRF-Residual * ! Also Mag.Sensor offset from GSP can be corrected * ! (for non-polar region data) * ! This program assumes no survey flight over midnight. * !************************************************************* implicit none real,parameter :: xdeg = 111111.1 ! length (m) of 1 deg. meridian real(8),parameter :: pai = (3.14159265359/180.) integer :: lwkdir, lrtrim real :: fid(300), fx(300), fy(300), fz(300) real(8) :: tmag(300), wtm(300), secs(300) character(80) :: wdr, flog, finf, fnam, fout character(160) :: line; character(80) :: str; character(57) :: s character(20) :: sdtm; character(8) :: name integer :: ldr, ngn, intv, kam, nps, mps, is, it, iymd, kst, kintv, IER integer :: i, nh, nm, iq, nsat, ifd, kc real :: forwd, strbd, dwnwd, dir, xnc, yec, zuc, ydeg, year, sec real :: falt, falt1, faltn, f, tigrf, tres real(8) :: tms, tme, tmsr, wtmr real(8) :: sec1, secn, secst, seced real(8) :: dido, dkdo, dido1, dkdo1, didon, dkdon, didop call premsg('XLDAM : Extract Line data from daq2asc output') call premsg(' with real-time Pos. data interporation') call premsg(' and calculation of IGRF-Residual') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr; finf = wdr; fnam = wdr; fout = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call gparma('Enter Lines.Inf. filename ==> ', 81-ldr, finf(ldr:80)) open(10,file=finf,status='old') call gparma('Enter Output filename ==> ', 81-ldr, fout(ldr:80)) open(20,file=fout,status='new') call gparmf(' Year of survey : ', year) call gparmi(' Generation number of IGRF : ', ngn) call gparmi(' Thin-out Ratio N (crop 1 out of N) ==> ', intv) if (intv <= 0) intv = 1 !*** call premsg(' Correction for Mag.Sensor offset from GPS antenna') call premsg(' (To disable correction, Enter blank line or "0 0 0")') call parmin(' [Forward, Starboard, Downward] (in m) ==> ', 80, str) kc = 1; forwd = -999.; strbd = -999.; dwnwd = -999. didop = 0.; ydeg = xdeg read(str,*,iostat=IER) forwd, strbd, dwnwd if ((forwd == -999.) .and. (strbd == -999.) .and. (dwnwd == -999.)) then forwd = 0.; strbd = 0.; dwnwd = 0. endif if ((forwd > 50.) .or. (forwd < -50.) .or. & (strbd > 50.) .or. (strbd < -50.) .or. & (dwnwd > 50.) .or. (dwnwd < -50.)) & call abendm('Abnormal offset data') if ((forwd == 0.) .and. (strbd == 0.) .and. (dwnwd == 0.)) kc = 0; !*** call clspin() write(6,'(a)') '----------------------------------------' write(6,'(a)') 'XLDAM : Extract Line data from daq2asc output' write(6,'(8x,a)') ' with real-time Pos. data interporation' write(6,'(8x,a)') ' and calculation of IGRF-Residual' write(6,'(a,a)') ' Lines.Inf. file : ', finf(1:lrtrim(finf)) write(6,'(a,a)') ' Output filename : ', fout(1:lrtrim(fout)) write(6,'(a,f7.2)') ' Year of survey : ', year write(6,'(a,i2)') ' IGRF generation : ', ngn write(6,'(a,i3)') ' Data thin-out ratio :', intv !*** write(6,'(a)') ' Mag.Sensor offset from GPS antenna (m)' write(6,'(3(6x,a,f6.2))') & 'Forward:', forwd, 'Starboard:', strbd, 'Downward:', dwnwd !*** call gigrf(ngn, year) open(9,status='scratch') kam = 0 do read(10,'(a)',iostat=IER) str if (IER < 0) exit if (str(1:1) == '=') then is = 2 do while ((is <= 80) .and. (str(is:is) == ' ')); is = is + 1; enddo if (is > 80) call abendm('Lines.Inf. data illegal') it = is + 1 do while ((it <= 80) .and. (str(it:it) /= ' ')); it = it + 1; enddo if (it > 80) call abendm('Lines.Inf. data illegal') fnam(ldr:80) = str(is:it-1) if (kam == 1) close(1); kam = 1 write(6,'(2a)') ' Input AM filename : ', fnam(1:lrtrim(fnam)) write(9,'(2a)') ' Input AM filename : ', fnam(1:lrtrim(fnam)) open(1,file=fnam,status='old') do read(1,'(a)') line if (line(1:14) == '//PC-Time Data') then write(6,'(2x,a)') line(1:45) write(9,'(2x,a)') line(1:45) endif if (line(1:10) == '/DateTime:') read(line,'(10x,i9)') iymd if (line(1:10) == '/ FID ') exit enddo do ! find first valid POS read(1,'(a)',iostat=IER) line if (IER < 0) call abendm('EOF reached on AM input') read(line,'(f7.1,1x,i2,1x,i2,1x,f5.2,10x,f10.3,3f8.3,40x,a57)') & fid(1), nh, nm, sec, tmag(1),fx(1),fy(1),fz(1),s wtm(1) = nh*10000 + nm*100 + sec secs(1) = dble(nh*3600 + nm*60) + sec; if (index(s,'*') == 0) then read(s,*) sec1,dido1,dkdo1,falt1,iq,nsat if (iq >= 1) exit endif enddo do i=2,300 ! find next valid POS read(1,'(a)',iostat=IER) line if (IER < 0) call abendm('EOF reached on AM input') read(line,'(f7.1,1x,i2,1x,i2,1x,f5.2,10x,f10.3,3f8.3,40x,a57)') & fid(i), nh, nm, sec, tmag(i),fx(i),fy(i),fz(i),s wtm(i) = nh*10000 + nm*100 + sec secs(i) = dble(nh*3600 + nm*60) + sec; if (index(s,'*') == 0) then read(s,*) secn,didon,dkdon,faltn,iq,nsat if (secn <= sec1) call abendm('AM data UTC illegal sequence') if (iq >= 1) then nps = i; exit endif endif enddo if (iq <= 0) then write(str(1:8),'(f8.1)') fid(1) call premsg('--- GPS data lost after Fid.=' // str(1:8)) endif sec = secs(nps) - secn write(6,'(a,f8.2,a,f7.1)') & ' AM time-offset: ', sec, ' sec at Fid.=', fid(nps) write(9,'(a,f8.2,a,f7.1)') & ' AM time-offset: ', sec, ' sec at Fid.=', fid(nps) cycle endif if (kam == 0) call abendm('input file not selected') is = 1 do while ((is <= 80) .and. (str(is:is) == ' ')); is = is + 1; enddo if (is > 80) call abendm('Lines.Inf. data illegal') it = is + 1 do while ((it <= 80) .and. (str(it:it) /= ' ')); it = it + 1; enddo if (it > 80) call abendm('Lines.Inf. data illegal') name = str(is:it-1) if (kc == 0) then read(str(it:80),*) tms, tme else !*** read(str(it:80),*) tms, tme, dir dir = dir * pai xnc = forwd*cos(dir) - strbd*sin(dir) yec = forwd*sin(dir) + strbd*cos(dir) zuc = -dwnwd !*** endif write(str,'(a8,1x,i9,1x,f9.2,1x,f9.2)') name, iymd, tms, tme call premsg(' Line: ' // str(1:38)) write(20,'(a1,a)') '&', str(1:38) if ((tms < 0.) .or. (tme < 0.) .or. (tme <= tms)) & call abendm('illegal TIME spec') nm = int(tms/100.); sec = sngl(tms - nm*100.) nh = nm/100; nm = nm - nh*100 if ((nh > 23) .or. (nm > 59) .or. (sec >= 60.)) & call abendm('illegal TIME spec') secst = dble(nh*3600 + nm*60) + sec nm = int(tme/100.); sec = sngl(tme - nm*100.) nh = nm/100; nm = nm - nh*100 if ((nh > 23) .or. (nm > 59) .or. (sec >= 60.)) & call abendm('illegal TIME spec') seced = dble(nh*3600 + nm*60) + sec kst = 0; kintv = 1 do do mps=1,(nps-1) if (secs(mps) > seced) exit if (secs(mps) < secst) cycle kintv = kintv - 1 if (kintv == 0) then f = (secs(mps)-sec1) / (secn-sec1) dido = dido1 + (didon-dido1)*f dkdo = dkdo1 + (dkdon-dkdo1)*f falt = falt1 + (faltn-falt1)*f call igrfc(sngl(dido), sngl(dkdo), falt, tigrf) tres = real(tmag(mps) - tigrf) if (tres < -9999.98) tres = -9999.99 if (tres > +9999.98) tres = +9999.99 ifd = nint(fid(mps) * 10.) !*** [Correction for Mag.Sensor offset from GPS, with 0.02m accuracy ! and assumption of sphere earth, for latitude range of 80S-80N deg.] if ((dido < -80.) .or. (dido > 80.)) & call abendm('Polar region not supported') if (abs(dido-didop) > 0.175) then didop = dido; ydeg = xdeg * cos(dido * pai) endif dido = dido + xnc/xdeg; dkdo = dkdo + yec/ydeg; falt = falt + zuc !*** write(20,'(i8,1x,i8,1x,f9.2,1x,i2,1x,f11.7,1x,f12.7,1x,f7.2,$)') & ifd, iymd, wtm(mps), 7, dido, dkdo, falt write(20,'(2(1x,f8.2),3(1x,f7.3),1x,f9.2)') & tmag(mps),tres, fx(mps),fy(mps),fz(mps), secs(mps) wtmr = wtm(mps); kintv = intv if (kst == 0) then tmsr = wtm(mps); kst = 1 endif endif enddo if (secs(mps) > seced) exit fid(1) = fid(nps); secs(1) = secs(nps); tmag(1) = tmag(nps) fx(1) = fx(nps); fy(1) = fy(nps); fz(1) = fz(nps); wtm(1) = wtm(nps) sec1 = secn; dido1 = didon; dkdo1 = dkdon; falt1 = faltn do i=2,300 ! find next valid POS read(1,'(a)',iostat=IER) line if (IER < 0) call abendm('EOF reached on AM input') read(line,'(f7.1,1x,i2,1x,i2,1x,f5.2,10x,f10.3,3f8.3,40x,a57)') & fid(i), nh, nm, sec, tmag(i),fx(i),fy(i),fz(i),s wtm(i) = nh*10000 + nm*100 + sec secs(i) = dble(nh*3600 + nm*60) + sec; if (index(s,'*') == 0) then read(s,*) secn,didon,dkdon,faltn,iq,nsat if (secn <= sec1) call abendm('AM data UTC illegal sequence') if (iq >= 1) then nps = i; exit endif endif enddo if (iq <= 0) then write(str(1:8),'(f8.1)') fid(1) call premsg('--- GPS data lost after Fid.=' // str(1:8)) if (kst == 1) call abendm('unable to recover position') do do ! find first valid POS read(1,'(a)',iostat=IER) line if (IER < 0) call abendm('EOF reached on AM input') read(line,'(f7.1,1x,i2,1x,i2,1x,f5.2,10x,f10.3,3f8.3,40x,a57)') & fid(1), nh, nm, sec, tmag(1),fx(1),fy(1),fz(1),s wtm(1) = nh*10000 + nm*100 + sec secs(1) = dble(nh*3600 + nm*60) + sec if (index(s,'*') == 0) then read(s,*) secn,didon,dkdon,faltn,iq,nsat if (iq >= 1) exit endif enddo do i=2,300 ! find next valid POS read(1,'(a)',iostat=IER) line if (IER < 0) call abendm('EOF reached on AM input') read(line,'(f7.1,1x,i2,1x,i2,1x,f5.2,10x,f10.3,3f8.3,40x,a44)') & fid(i), nh, nm, sec, tmag(i),fx(i),fy(i),fz(i),s wtm(i) = nh*10000 + nm*100 + sec secs(i) = nh*3600 + nm*60 + sec if (index(s,'*') == 0) then read(s,*) secn,didon,dkdon,faltn,iq,nsat if (secn <= sec1) call abendm('AM data UTC illegal sequence') if (iq >= 1) then nps = i; exit endif endif enddo if (iq <= 0) then write(str(1:8),'(f8.1)') fid(1) call premsg('--- GPS data lost after Fid.=' // str(1:8)) endif enddo endif enddo if (kst == 0) call abendm('No data Extracted') write(str,'(a17,f9.2,1x,f9.2,a10,i4)') & ' Extracted : ', tmsr, wtmr, ' intv =', intv write(6,'(a)') str(1:50) write(9,'(a)') str(1:50) enddo close(1); close(10); close(20) endfile(9); rewind(9) if (flog(ldr:ldr) /= ' ') then open(99,file=flog,access='append') write(99,'(//a)') 'XLDAM : Extract Line data from daq2asc output' write(99,'(8x,a)') ' with real-time Pos. data interporation' write(99,'(8x,a)') ' and calculation of IGRF-Residual' write(99,'(a,a)') ' Lines.Inf. file : ', finf(1:lrtrim(finf)) write(99,'(a,a)') ' Output filename : ', fout(1:lrtrim(fout)) write(99,'(a,f7.2)') ' Year of survey : ', year write(99,'(a,i2)') ' IGRF generation : ', ngn write(99,'(a,i3)') ' Data thin-out ratio :', intv if (kc == 0) then write(99,'(a)') ' Not corrected for Mag.Sensor offset from GPS' else write(99,'(a)') ' Corrected for Mag.Sensor offset (m) by' write(99,'(3(6x,a,f6.2))') & 'Forward:', forwd, 'Starboard:', strbd, 'Downward:', dwnwd endif do read(9,'(a)',iostat=IER) line if (IER < 0) exit write(99,'(a)') line(1:lrtrim(line)) enddo call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif close(9) stop end