!************************************************************* ! XLDPN : Extract Line data from daq2asc output * ! with PostNav 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 character(80) :: wdr, flog, finf, fnam, fpos, fout character(160) :: line; character(80) :: str; character(60) :: ss character(20) :: sdtm; character(8) :: name; character(5) :: tzone integer :: ldr, ngn, intv, kam, is, it, kpn, iymd, kst, kintv, IER integer :: k, ih, im, ifd, kc, kt, mintz real :: forwd, strbd, dwnwd, dir, xnc, yec, zuc, ydeg real :: fid, fx, fy, fz, year, sec real :: falt, falta, faltb, f, tigrf, tres real(8) :: wtm, tms, tme, tmsr, tmer real :: sect, seca, secb, secst, seced real(8) :: tmag, dido, didoa, didob, dkdo, dkdoa, dkdob, didop call premsg('XLDPN : Extract Line data from daq2asc output') call premsg(' with PostNav Pos. data interporation') call premsg(' and calculation of IGRF-Residual') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr; finf = wdr; fnam = wdr; fpos = 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 premsg(' Time Spec. of Position data ?') call gparmi(' Localtime sama as Obs.Mag.data [0] or UTC[1] ==> ', kt) if (kt == 0) then mintz = 0 else if (kt == 1) then call premsg(' UTC time will be converted into Localtime.') call gparma(' Localtime-Zone ? "[+/-]HHMM" (default: +0900) ==> ', & 5, tzone) if (tzone == ' ') tzone = '+0900' if ((tzone(1:1) /= '+') .and. (tzone(1:1) /= '-')) & call abendm('invalid timezone specification') read(tzone(2:5),'(2i2)') ih, im if ((ih<0) .or. (ih>23) .or. (im<0) .or. (im>59)) & call abendm('invalid timezone specification') mintz = ih*60 + im if (tzone(1:1) == '-') mintz = -mintz else call abendm('invalid specification') endif 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)') 'XLDPN : Extract Line data from daq2asc output' write(6,'(8x,a)') ' with PostNav Pos. data interporation' write(6,'(8x,a)') ' and calculation of IGRF-Residual' write(6,'(2a)') ' Lines.Inf. file : ', finf(1:lrtrim(finf)) if (kt == 0) then write(6,'(a)') ' Input Pos. data Time : Localtime (same as Obs.Mag.)' else write(6,'(a)') ' Input Pos. data Time : UTC' write(6,'(a,a)') ' will be converted into Localtime of Zone=', tzone endif write(6,'(2a)') ' 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; kpn = 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) is = it + 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') fpos(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 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)') & fid, ih, im, sec, tmag,fx,fy,fz wtm = ih*10000 + im*100 + sec sect = dble(ih*3600 + im*60) + sec; if (kpn == 1) close(2); kpn = 1 write(6,'(2a)') ' PNAV data filename : ', fpos(1:lrtrim(fpos)) write(9,'(2a)') ' PNAV data filename : ', fpos(1:lrtrim(fpos)) open(2,file=fpos,status='old') read(2,'(a)',iostat=IER) ss ! get first PNAV data if (IER < 0) call abendm('EOF reached on PNAV input') k = index(ss, ':') if ((k == 0) .or. (ss(k+3:k+3) /= ':') .or. (ss(k+6:k+6) /= '.')) & call abendm('PNAV data illegal format') ss(k:k) = ' '; ss(k+3:k+3) = ' ' read(ss,*) ih, im, sec, didoa, dkdoa, falta if ((ih < 0) .or. (im < 0) .or. (sec < 0.) .or. & (ih >= 24) .or. (im >= 60) .or. (sec >= 60.)) & call abendm('PNAV data illegal format') seca = dble((ih*60 + im + mintz) * 60) + sec if (seca >= 86400.) seca = seca - 86400. if (seca < 0.) seca = seca + 86400. read(2,'(a)',iostat=IER) ss ! get next PNAV data if (IER < 0) call abendm('EOF reached on PNAV input') k = index(ss, ':') if ((k == 0) .or. (ss(k+3:k+3) /= ':') .or. (ss(k+6:k+6) /= '.')) & call abendm('PNAV data illegal format') ss(k:k) = ' '; ss(k+3:k+3) = ' ' read(ss,*) ih, im, sec, didob, dkdob, faltb secb = dble((ih*60 + im + mintz) * 60) + sec if (secb >= 86400.) secb = secb - 86400. if (secb < 0.) secb = secb + 86400. if (secb <= seca) call abendm('PNAV data [sec] illegal sequence') 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 read(str(it:80),*) tms, tme 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') im = int(tms/100.); sec = sngl(tms - im*100.) ih = im/100; im = im - ih*100 if ((ih > 23) .or. (im > 59) .or. (sec >= 60.)) & call abendm('illegal TIME spec') secst = dble(ih*3600 + im*60) + sec im = int(tme/100.); sec = sngl(tme - im*100.) ih = im/100; im = im - ih*100 if ((ih > 23) .or. (im > 59) .or. (sec >= 60.)) & call abendm('illegal TIME spec') seced = dble(ih*3600 + im*60) + sec kst = 0; kintv = 1 do while (sect < secst) 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)') & fid, ih, im, sec, tmag,fx,fy,fz wtm = ih*10000 + im*100 + sec sect = dble(ih*3600 + im*60) + sec enddo if (sect > seced) call abendm('Nodata Extracted') do while (wtm <= tme) if (seca > sect) then write(ss,'(a10,f9.1,a26,f9.1)') & 'PNAV (sec=', seca, ') passed AMasc time [sec=]', sect call abendm(ss) endif do while (secb < sect) seca = secb; didoa = didob; dkdoa = dkdob; falta = faltb read(2,'(a)',iostat=IER) ss ! get next PNAV data if (IER < 0) call abendm('EOF reached on PNAV input') k = index(ss, ':') if ((k == 0) .or. (ss(k+3:k+3) /= ':') .or. (ss(k+6:k+6) /= '.')) & call abendm('PNAV data illegal format') ss(k:k) = ' '; ss(k+3:k+3) = ' ' read(ss,*) ih, im, sec, didob, dkdob, faltb secb = dble((ih*60 + im + mintz) * 60) + sec if (secb >= 86400.) secb = secb - 86400. if (secb < 0.) secb = secb + 86400. if (secb <= seca) then write(ss,'(a10,f9.1,a21,f9.1)') & 'PNAV (sec=', secb, ') illegal seq. after ', seca call abendm(ss) endif enddo kintv = kintv - 1 if (kintv == 0) then f = (sect-seca) / (secb-seca) dido = didoa + (didob-didoa)*f dkdo = dkdoa + (dkdob-dkdoa)*f falt = falta + (faltb-falta)*f call igrfc(sngl(dido), sngl(dkdo), falt, tigrf) tres = real(tmag - tigrf) if (tres < -9999.98) tres = -9999.99 if (tres > +9999.98) tres = +9999.99 ifd = nint(fid * 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,'(2(i8,1x),f9.2,1x,i2,1x,f11.7,1x,f12.7,1x,f7.2,$)') & ifd, iymd, wtm, 3, dido, dkdo, falt write(20,'(2(1x,f8.2),3(1x,f7.3),1x,f9.2)') & tmag, tres, fx, fy, fz, sect tmer = wtm; kintv = intv if (kst == 0) then tmsr = wtm; kst = 1 endif endif 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)') & fid, ih, im, sec, tmag,fx,fy,fz wtm = ih*10000 + im*100 + sec sect = dble(ih*3600 + im*60) + sec enddo if (kst == 0) call abendm('No data Extracted') write(str,'(a17,f9.2,1x,f9.2,a10,i4)') & ' Extracted : ', tmsr, tmer, ' intv =', intv write(6,'(a)') str(1:50) write(9,'(a)') str(1:50) enddo close(1); close(2); close(10); close(20) endfile(9); rewind(9) if (flog(ldr:ldr) /= ' ') then open(99,file=flog,access='append') write(99,'(//a)') 'XLDPN : Extract Line data from daq2asc output' write(99,'(8x,a)') ' with PostNav Pos. data interporation' write(99,'(8x,a)') ' and calculation of IGRF-Residual' write(99,'(2a)') ' Lines.Inf. file : ', finf(1:lrtrim(finf)) if (kt == 0) then write(99,'(a)') ' Input Pos. data Time : Localtime (same as Obs.Mag.)' else write(99,'(a)') ' Input Pos. data Time : UTC' write(99,'(a,a)') ' were converted into Localtime of Zone=', tzone endif write(99,'(2a)') ' 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