program samla
! Originally written by Gunstein Dalane at Kartverket, Norway
! Edited:
! 2017-03-07: Forced name sorting in DIR command. Adding comments. Changed format for reading 
!             coordinates to 3f15.5 (Bernese ver 5.2)./Lotti
! 2017-03-09: Computation of average and differences to averge, converted to meter. Year converted to decimal year. /Lotti
! 2017-03-17: Modified output format for Tsview and added file for Hector/Lotti
! 2017-03-21: Changed names of the mb-files for TSview /Lotti
! 2017-04-24: Produce enu-files instead of NEU-files. (Contineous time scale with MDJ, differences in mm)/Lotti
! 2017-04-26: Sort the years 1900 before 2000. Wrongly sorted before as we just use 2 digit years in the file
!             names. /Lotti
! 2018-01-31: Changed back to compute the decimal years by dividing doy by 365 or 366, works better with TSview.
!             Changed also the formula for computing MJD. /Lotti

implicit none
integer*4 fil_l,i,j,pkt_l,ii,iii,bgr,bmin,lgr,lmin,istat,ant,stat,intnum,inr,ir
integer*4 antal_crd(10000),ndays,n_days,iiii,iii2, iii3, iii4,yy,yyyy, ddd   
integer*4 n_leap, mjd 
character fil*20,string*80,pkt*4,flag*1,stasjonsfil*8,pnr(1000)*4,stasjonsfil_enu*12
character stasjonsfil_n*16,stasjonsfil_e*16,stasjonsfil_u*16
real*8 x,y,z,wgs84_a,wgs84_e2,wgs84_f,br,rl,h,bsek,lsek, dlat, dlong,dn,de,dh,yyyydec 
real*8  x_m(1000), y_m(1000), z_m (1000), lat_m(1000),long_m(1000),h_m(1000),pi
real*8 mjddec
logical*4 open
data wgs84_a,wgs84_e2,wgs84_f/6378137.0D0,0.00669437999013D0,0.00335281066474/
data pi/3.1415926535897932d0/
call overskrift
call system("del *.dat")
call system("del *.dat1")
call system("del *.dat2")
call system("del *.dat3")
call system("del *.enu")
call system("dir /B /O N *.CRD > LSFIL")

! Sort the LSFIL so yy >89 (1990-99) will come before years < 90 (2000-89)
open(7,file = 'LSFIL',status = 'old',access = 'sequential',iostat = istat)
open (8, file = 'LSFIL_SORT', access = 'sequential',iostat = istat)
fil_l = 2
do while(fil_l >= 2)
  read(7,'(a)',iostat = istat,end = 701)fil
  fil_l = len_trim(fil)
  read(fil(4:5),*) yy  
  if (yy >89) then 
     write (8,'(a)') fil
  end if
end do
701 continue  
rewind (7)
fil_l = 2
do while(fil_l >= 2)
  read(7,'(a)',iostat = istat,end = 702)fil
  fil_l = len_trim(fil)
  read(fil(4:5),*) yy  
  if (yy <90) then 
     write (8,'(a)') fil
  end if
end do
702 continue
close (7)
close (8)

! Read through the list of coordinate files, populate the array pnr with all station names (4ID) and count nr of stn
! as well as average for each station
open(8,file = 'LSFIL_SORT',status = 'old',access = 'sequential',iostat = istat)
fil_l = 2
i = 0
ant = 0
do while(fil_l >= 2)
  read(8,'(a)',iostat = istat,end = 901)fil
  fil_l = len_trim(fil)
  open(10,file = fil,status = 'old',access = 'sequential',iostat = istat,err = 901)  !For each coordinate file
  do j = 1,6
    read(10,'(a)',end = 801)string
  end do
  pkt_l = 2
  do while(pkt_l >=  2)
    read(10,'(t1,i3,t6,a,t22,3f15.5,t71,a)',iostat = istat,end = 801)ii,pkt,x,y,z,flag
    if(flag == 'A' .or. flag == 'W')then
      call tolower(pkt)
      inr = intnum(pkt,ant,pnr)
      if(inr.ne.0)go to 998     ! Station alreday exist in pnr
      i = i + 1
      ant = i
      pnr(i) = pkt
      x_m(i)=0
      y_m(i)=0
      z_m(i)=0
      antal_crd(i)=0
      inr=i
!998   continue
998   x_m(inr)=x+x_m(inr)   
      y_m(inr)=y+y_m(inr)
	  z_m(inr)=z+z_m(inr)
      antal_crd(inr)=antal_crd(inr)+1
    end if
  end do
801 continue
end do
901 continue



open(99,file = 'Average.txt',access = 'sequential',iostat = istat)

do j = 1,ant                                  ! Compute average for each station and convert to lat and long
    x_m(j)=x_m(j)/antal_crd(j)  
    y_m(j)=y_m(j)/antal_crd(j)
    z_m(j)=z_m(j)/antal_crd(j)
    lat_m(j)=0.0d0
    call bxtran(wgs84_a,wgs84_f,lat_m(j),long_m(j),h_m(j),x_m(j),y_m(j),z_m(j),ir)
	write(99,'(t1,a,t7,3(f20.4),2(f20.11),f10.4)',iostat = stat) pnr(j),x_m(j),y_m(j),z_m(j),lat_m(j)*180/pi,long_m(j)*180/pi,h_m(j)
end do
close (99)
rewind (8)
close (10)

! Read through the list of coordinate files again and get coordinates, convert them to lat,long, ellh and differencies
fil_l = 2
do while(fil_l >= 2)
  read(8,'(a)',iostat = istat,end = 900)fil
  fil_l = len_trim(fil)
  open(10,file = fil,status = 'old',access = 'sequential',iostat = istat,err = 900)
  do j = 1,6
    read(10,'(a)',end = 800)string
  end do
  pkt_l = 2
  do while(pkt_l >=  2)
    read(10,'(t1,i3,t6,a,t22,3f15.5,t71,a)',iostat = istat,end = 800)ii,pkt,x,y,z,flag
    if(flag == 'A' .or. flag == 'W')then
      call tolower(pkt)
      br = 0.0d0
      pkt_l = len_trim(pkt)
      call bxtran(wgs84_a,wgs84_f,br,rl,h,x,y,z,ir)
      call r2gms(br,bgr,bmin,bsek)
      call r2gms(rl,lgr,lmin,lsek)
      inr = intnum(pkt,ant,pnr)

      dlat=br-lat_m(inr)                                      !Compute differences to average
      dlong=rl-long_m(inr)
      dh=h-h_m(inr)

      call lld2m(wgs84_a,wgs84_e2,br,rl,dlat,dlong,dn,de)     !Convert differencies

      read(fil(4:5),*) yy                                     ! Convert year to integer
      read(fil(6:8),*) ddd                                    ! Convert day to integer
	  
	  if(yy>90) then
	     yyyy=1900.0+yy
      else
         yyyy=2000.0+yy
      end if

      if(yyyy/4.0 - INT(yyyy/4.0) < 0.1) then                 ! Check if leap year
         n_days=366
      else 
         n_days=365
      end if		 
      yyyydec = REAL(yyyy) + (REAL(ddd)-0.5)/REAL(n_days)     ! Convert to decimal year, maybe not becessary to be put REAL everywhere but safest 

! Compute MJD from yyyy and ddd
      mjddec=REAL(mjd(yyyy,ddd)) + 0.5                ! Add 0.5 day to get time stamp at noon (12:00)

! Compute decimal year from MJD
! Julian epoch J = 2000.0 + (MDJ - 51544.5)/365.25     to be sure to always be within the correct day .5 can be changed to .4
!      yyyydec=(mjddec-51544.4)/365.25 + 2000            ! same difference between all days , NOT USED for TSview !!!
	  
      iii = inr + 100  ! ii i= inr + 100              ! File handle for stationwise file with lat, long, height
      inquire(unit = iii,opened = open)      
      stasjonsfil = pkt//'.dat'
      if(open .eqv. .false.)then
        open(iii,file = stasjonsfil,access = 'sequential',iostat = istat)
      end if
      write(iii,'(t1,a,t7,2(2i5,f11.5),f11.4)',iostat = stat)fil(4:5)//'/'//fil(6:9),bgr,bmin,bsek,lgr,lmin,lsek,h

      iiii = inr + 1000                 ! File handle for stationwise file with dE, dN, dU
      inquire(unit = iiii,opened = open)      
      stasjonsfil_enu = pkt//'_'//fil(1:3)//'.enu'
      if(open .eqv. .false.)then
        open(iiii,file = stasjonsfil_enu,access = 'sequential',iostat = istat)
      end if
      write(iiii,'(t1,f7.1,t12,f10.2,f10.2,f10.2)',iostat = stat)mjddec,de*1000, dn*1000, dh*1000

      iii2 = inr + 2000                 ! File handle for stationwise file with dN-diff for TsView
      inquire(unit = iii2,opened = open)      
      stasjonsfil_n = 'mb_'//pkt//'_GPS.dat1'
      if(open .eqv. .false.)then
        open(iii2,file = stasjonsfil_n,access = 'sequential',iostat = istat)
		write(iii2,'(a)',iostat = stat) 'Bernese analysis '//fil(1:3)//'-solutions'
		write(iii2,'(t1,a,a,t29,f13.3,a)',iostat = stat) pkt,'_GPS to N Solution  1 + ', lat_m(inr)*20000000/pi, ' m'
		write(iii2,'(a)', iostat = stat) ''
      end if
      write(iii2,'(t1,f10.5,t13,f12.5,a)',iostat = stat) yyyydec,dn,'  0.00100'

      iii3 = inr + 3000                 ! File handle for stationwise file with dE-diff for TsView
      inquire(unit = iii3,opened = open)      
      stasjonsfil_e = 'mb_'//pkt//'_GPS.dat2'
      if(open .eqv. .false.)then
        open(iii3,file = stasjonsfil_e,access = 'sequential',iostat = istat)
		write(iii3,'(a)',iostat = stat) 'Bernese analysis '//fil(1:3)//'-solutions'
		write(iii3,'(t1,a,a,t29,f13.3,a)',iostat = stat) pkt,'_GPS to E Solution  1 + ', long_m(inr)*cos(lat_m(inr))*20000000/pi, ' m'
		write(iii3,'(a)', iostat = stat) ''		
      end if
      write(iii3,'(t1,f10.5,t13,f12.5,a)',iostat = stat) yyyydec,de,'  0.00100'

      iii4 = inr + 4000                 ! File handle for stationwise file with dU-diff for TsView
      inquire(unit = iii4,opened = open)      
      stasjonsfil_u = 'mb_'//pkt//'_GPS.dat3'
      if(open .eqv. .false.)then
        open(iii4,file = stasjonsfil_u,access = 'sequential',iostat = istat)
		write(iii4,'(a)',iostat = stat) 'Bernese analysis '//fil(1:3)//'-solutions'
		write(iii4,'(t1,a,a,t29,f13.3,a)',iostat = stat) pkt,'_GPS to U Solution  1 + ', h_m(inr), ' m'
		write(iii4,'(a)', iostat = stat) ''		
		end if
      write(iii4,'(t1,f10.5,t13,f12.5,a)',iostat = stat) yyyydec,dh,'  0.00300'
	  
    end if   
  end do
800 continue
end do
900 continue

close (8)
!call system("del LSFIL")
end

!-------------------------------------------------------------------------------
!  SUBROUTINES and FUNCTIONS
!-------------------------------------------------------------------------------

subroutine tolower(text)
implicit none
character text*(*)
integer*2 i,j,l,lengde
l = len_trim(text)
do i = 1,l
  j = ichar(text(i:i))
  if(j >= 65 .and. j <= 93) text(i:i) = char(j + 32)
end do
return
end


SUBROUTINE BXTRAN(A,F,BR,RL,H,X,Y,Z,IR)
!        MELLOM BREDDE,LENGDE OG REKTANGUL[RE
!        BR=0:TIL BREDDE LENGDE
!        BRUKER BOWRINGS FORMEL:SURVEY REVIEW XXIII:181,pp323-327
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IR=0
E2=F*(2.D0-F)
IF(A.LT.6.D6.OR.A.GT.7.D6)IR=61
IF(F.LT.0.003D0.OR.F.GT.0.004D0)IR=62
IF(dabs(BR) < 1.D-20)GO TO 100
IF(BR.GT.1.58D0.OR.BR.LT.-1.58D0)IR=63
IF(RL.GT.3.15D0.OR.RL.LT.-3.15D0)IR=64
IF(IR.NE.0)RETURN
SIB=DSIN(BR)
RN=A/DSQRT(1.D0-E2*SIB**2)
P=(RN+H)*DCOS(BR)
X=P*DCOS(RL)
Y=P*DSIN(RL)
Z=(RN*(1.D0-E2)+H)*SIB
RETURN
100 P=DSQRT(X**2+Y**2)
IF(DABS(P).LT.0.0001)P=0.0001
TQ=Z/(P*(1.D0-F))
C2=1.D0/(1.D0+TQ**2)
CQ=DSQRT(C2)
!      Q=DATAN(TQ)
!      CQ=DCOS(Q)
SQ=CQ*TQ
CQ3=C2*CQ
SQ3=(1.D0-C2)*SQ
IF(DABS(X).LT.0.0001)X=0.0001
!      RL=DATAN(Y/X)
RL=DATAN2(Y,X)
TA=(Z+E2*A*SQ3/(1.D0-F))/(P-E2*A*CQ3)
BR=DATAN(TA)
T2=TA**2
C2I=1.D0+T2
S2=T2/C2I
RN=A/DSQRT(1.D0-E2*S2)
H=P*DSQRT(C2I)-RN
RETURN
END

subroutine r2gms(rad,gr,min,sek)
integer*4 gr,min
real*8 rad,sek,ro,dgr,dmin
ro = 180.0d0/3.1415926535897932d0
dgr = rad*ro
gr = dgr
dmin = (dgr-dfloat(gr))*60.0d0
min = dmin
sek = (dmin-dfloat(min))*60.0d0
return
end 

! Subroutine to convert differences in lat and longitude to meter
! Formler för meridiankrökningsradie och tvärkrökningsradie ur HMK_Ge:S, Bilaga B
! lld2m (=lat long differences to meter), by Lotti Jivall
subroutine lld2m(a,e2,lat,long,dlat,dlong,dn,de)
real*8 a,e2,lat,long,dlat,dlong,dn,de,W,M,N
W = SQRT(1-e2*sin(lat)*sin(lat))
M = a*(1-e2)/W/W/W
N = a/W
dn=M*dlat
de=N*cos(lat)*dlong
return
end 


subroutine overskrift
! call system("clear")
write(*,'(17(t5,a,/))')&
'*******************************************************************',&
'*                                                                 *',&
'*                 P R O G R A M   S A M L A                       *',&
'*                                                                 *',&
'*           is a modification of the program SAMLE                *',& 
'*           written by Gunstein Dalane, Kartverket, Norway        *',&
'*            modified by Lotti Jivall, Lantmateriet, Sweden       *',&
'*                                                                 *',&
'*          T H I S  P R O G R A M   E X T R A C T S               *',&
'*             X Y X - COORDINATES FROM *.CRD FILES                *',&
'*         and compiles (samler/samlar) them station by station    *',&
'*         as latitude, longtitude and ellipsoidal heights         *',&
'*      as well as three stationwise files with differences in     *',&
'*     north, east and up for input to TsView and one .enu-file    *',&
'*                     for input to Hector                         *',&
'*                                                                 *',&
'*******************************************************************' 
return
end

! ******* Function som finner et punktnummer's plassering i pulja. ******
integer*4 function intnum(gradnr,pant,pytrnr)
integer*4 pant,i
character gradnr*(*),pytrnr(1:pant)*(*)
intnum=0
i = 1
do while(i.le.pant)
  if(gradnr.eq.pytrnr(i))then
    intnum=i
    i=pant+1
  end if
  i = i+1
end do
return
end

!! WRONG ******* Function to convert from year and day to MJD. ****** Not 100% correct!!
!integer*4 function mjd(year,day)
!integer*4 year,day,n_leap
!
!! number of leapdays since MDJ=0 (1858-11-17) 
!! note that it is first the year after a leapyear that the extra day has to be counted 
!n_leap=INT((year-1857)/4+0.00001) - INT((year-1801)/100+0.00001) + INT((year-1601)/400+0.01)
!
!mjd=(year-1858)*365 + n_leap + day  - 321
!return 
!end

! Function to convert from year and day (doy) to MJD. ******
! Formulas from Hofmann&Wellenhof adjusted to MDJ and doy.
! MJD=JD - 2400000.5
! JD=INT(365.25*y) + INT(30.6001(m+1)) D UT/24 + 1720981.5
! If M<=2: y=Y-1, m=M+12
! If M>2: y=Y, m=M
! -> MJD=INT(365.25(Y-1)) + INT(30.60001*14) + doy - 679019



integer*4 function mjd(year,day)
integer*4 year,day

mjd=INT(365.25*(year-1)) + INT(30.60001*14) + day - 679019 
return 
end
