      program setinput
C................................................................Aug 2019
C Insert TECgps to IRI-Plas input file
C Produce fitting coefficients for NmF2(TEC) and hmF2(TEC)
C
      dimension x(24),y(24),ynm(24)
      REAL tec(0:23),fm(0:23),hm(0:23),fp(0:23),hp(0:23),xnm(0:23)
C	INTEGER idatf(0:23),idath(0:23),idatt(0:23)
      character*2 AYR,AMN,ADY,BYR,BMN,BDY,ST
      character*68 infilef,infileh,infilet,outfile,outfilef,outfileh
	character*4 stat,jjjj
	character*3 jjj
C
      jjjj='0000'
	jjj='0.0'
	j1=0
	j2=-1
C
	outfile='statYRMNDY.inp'
      outfilef='polyymmddfp.txt'   ! Mathlab polyfit2
      outfileh='polyymmddhp.txt'   ! Mathlab polyfit2

	write(*,*) ' Enter AYR, AMN, ADY:'
	read(*,*) AYR, AMN, ADY
   30 FORMAT(3A2,1X,24(F4.1,1X))
   31 FORMAT(3A2,A68,F6.2)
   33 FORMAT(3A2,1X,24(F4.0,1X))
	outfile(5:6)=AYR
	outfile(7:8)=AMN
	outfile(9:10)=ADY
C
	outfilef(4:5)=AYR
	outfileh(4:5)=AYR
	outfilef(6:7)=AMN
      outfileh(6:7)=AMN
	outfilef(8:9)=ADY
      outfileh(8:9)=ADY
      OPEN(14,file=outfilef,access='APPEND')		! fit foF2 coefs
      OPEN(15,file=outfileh,access='APPEND')	 	! fit hmF2 coefs


	infilet='d:\web\graf\20YR\stYRt.txt'
	infilet(15:16)=AYR
	infilet(20:21)=AYR
	infilef=infilet
	infilef(22:22)='f'
	infileh=infilet
	infileh(22:22)='h'
    1 write(*,*) ' Enter ST:'
      read(*,*) ST
	write(*,*) ' Enter STAT or <0000> STOP:'
      read(*,*) STAT
	if (stat.eq.jjjj) GOTO 22  ! STOP
      write(*,*) ' Enter glat,glon:'
	read(*,*) glat,glon
C
	outfile(1:4)=STAT
	OPEN(10,file=outfile) ! .inp file for IRI-Plas
C
	infilef(18:19)=st
	infileh(18:19)=st
	infilet(18:19)=st
      OPEN(11,file=infilef,action='READ')
	OPEN(12,file=infileh,action='READ')
	OPEN(13,file=infilet,action='READ')
    2	read(13,30,end=40) BYR,BMN,BDY,(tec(k),k=0,23)
      read(11,30,end=40) BYR,BMN,BDY,(fm(k),k=0,23)
	read(12,33,end=40) BYR,BMN,BDY,(hm(k),k=0,23)
       IF ((AYR.eq.BYR).and.(AMN.eq.BMN).and.(ADY.eq.BDY)) then
	goto 40
	 ELSE
	goto 2
	 ENDIF
   40	close(unit=11)
      close(unit=12)
	close(unit=13)
C Prepare fitting foF2:
      DO i=0,23
	x(i+1)=tec(i)
	ynm(i+1)=1.24*fm(i)**2
	y(i+1)=ynm(i+1)  ! First NmF2(TEC)
	ENDDO

C Produce p1f,p2f,p3f for NmF2(TEC)
      call subpoly3c(x,y,p1f,p2f,p3f,rmse)
	write(*,*) ' NmF2 vs TEC'
C+      write(13,*) ' NmF2 vs TEC'	 
	write(*,*) stat,' p1= ', p1f, ' p2= ',p2f,' p3=',p3f,' rms=',rmse
	write(14,*) stat, p1f, p2f, p3f
C Output profile:
 	DO jj=0,23
	xx=tec(jj)
	xnm(jj)=xx*(p1f*xx+p2f)+p3f   ! NmF2
	if (xnm(jj).lt.4.0) xnm(jj)=ynm(jj+1)
	fp(jj)=sqrt(xnm(jj)/1.24)
	write(*,*) jj,xx,fp(jj)
C+	write(11,*) jj,xx,fp(jj)
	ENDDO

C Prepare fitting hmF2:
      DO i=0,23
	y(i+1)=hm(i)  ! Second hmF2(TEC)
	ENDDO

C Produce p1h,p2h,p3h for hmF2(TEC)
      call subpoly3c(x,y,p1h,p2h,p3h,rmse)
	write(*,*) ' hmF2 vs TEC'
C+      write(14,*) ' hmF2 vs TEC'	 
	write(*,*) stat,' p1= ', p1h, ' p2= ',p2h,' p3=',p3h,' rms=',rmse
	write(15,*) stat, p1h, p2h, p3h
C Output profile:
 	DO jj=0,23
	xx=tec(jj)
	hp(jj)=xx*(p1h*xx+p2h)+p3h   ! NmF2
	if ((hp(jj).gt.555.0).or.(hp(jj).lt.175.0)) hp(jj)=hm(jj+1)
	write(*,*) jj,xx,hp(jj)
C+	write(11,*) jj,xx,fp(jj)
	ENDDO

   32 FORMAT(3A2,1X,A4,1X,F4.1,3X,I1,2X,I2,1X,I2,1X,I1,2F6.1,3(3X,A3)
     +,2X,I4,1X,F4.1,4X,F6.2)
C Prepare .inp file first out from /graf/YEAR/stYRf.txt,/stYRh.txt,/stYRt.txt,
	DO iut=0,23
	rut=float(iut)+25.0
	ihmF2=nint(hm(iut))
	write(*,32) AYR,AMN,ADY,jjjj,rut,j1,j2,j2,j1,glat,glon,jjj,jjj
	+,jjj,ihmF2,fm(iut),tec(iut)
      write(10,32) AYR,AMN,ADY,jjjj,rut,j1,j2,j2,j1,glat,glon,jjj,jjj
	+,jjj,ihmF2,fm(iut),tec(iut)
	ENDDO
C Prepare .inp file Then proxy fpF2 and hmF2 (TEC=0.0)
	DO iut=0,23
	rut=float(iut)+25.0
	ihmF2=nint(hp(iut))
	zertec=0.
	write(*,32) AYR,AMN,ADY,jjjj,rut,j1,j2,j2,j1,glat,glon,jjj,jjj
	+,jjj,ihmF2,fp(iut),zertec
      write(10,32) AYR,AMN,ADY,jjjj,rut,j1,j2,j2,j1,glat,glon,jjj,jjj
	+,jjj,ihmF2,fp(iut),zertec
	ENDDO
   41	close(unit=10)
C
      write(*,*) STAT
	      GOTO 1
   22	      close(unit=14)
	      close(unit=15)
	pause ' '
	stop     
	END
C=========================================
      subroutine subpoly3c(x,y,p1,p2,p3,rmse)
	dimension x(24),y(24) 
	nut=24  ! diurnal UT set
C  
	s0=24.
	s1=0.
	s2=0.
	s3=0.
	s4=0.
	syx2=0.
	syx=0.
	sy=0.
      DO N=1,nut
	x1=x(N)
	x2=x1*x1
	x3=x1*x2
	x4=x2*x2
	y1=y(N)
	yx=y1*x1
	yx2=y1*x2
	s1=s1+x1
	s2=s2+x2
	s3=s3+x3
	s4=s4+x4
	syx2=syx2+yx2
	syx=syx+yx
	sy=sy+y1
      ENDDO 
C
      d1=(s0*s3-s1*s2)
      c1=(s0*s4-s2*s2)/d1
      c2=(s0*syx2-sy*s2)/d1
	d2=(s0*s2-s1*s1)
      c3=(s0*s3-s1*s2)/d2
      c4=(s0*syx-sy*s1)/d2

C
C p1=a; p2=b; p3=c
C
C polyfit2 coefficients
	p1=(c2-c4)/(c1-c3)
	p2=c2-p1*c1
	p3=(sy-p1*s2-p2*s1)/s0
C RMSE
      sum=0.
	DO N=1,nut
	del=y(N)-(x(N)*(p1*x(N)+p2)+p3)
	sum=sum+del**2
	ENDDO
	rmse=sqrt(sum/nut)
C
	RETURN
	END
C

