C+
C NAME:
C       program ashi_ratio
C PURPOSE:
C       To read the distances from the image axis and determine the distance in degrees relative to the number of pixels
C
C CALLING SEQUENCE:
C       ./ashi_ratio
C INPUTS:
C     N           integer  number of distance determinations = 4
C     distx(N,I)  real     X distance in pixels from the axis
C     disty(N,I)  real     Y distance in pixels from the axis
C     offset1     real     offset for determination 1    
C     offset2     real     offset for determination 2
C     offset3     real     offset for determination 3
C     offset4     real     offset for determination 4
C
C FUNCTIONS/SUBROUTINES:
C     read_ashi_dist
C
C MODIFICATION HISTORY:
C     November 2022, Bernard Jackson (UCSD)
C-    

        program ashi_ratio
      
          parameter(N =4)               ! Number of determinations
          parameter(INUM =25)           ! Number of x, y values
      
          real  distx(N,INUM),          ! X distance in pixels from the axis
     &          disty(N,INUM)           ! Y distance in pixels from the axis

          real  disttp1(INUM),
     &          disttp2(INUM),       
     &          disttp3(INUM),
     &          disttp4(INUM),
     &          degree1(INUM),
     &          degree2(INUM),
     &          degree3(INUM),
     &          degree4(INUM)

           character cfile1*18 /'ruler_points_1.txt'/,
     &               cfile2*18 /'ruler_points_2.txt'/,           
     &               cfile3*18 /'ruler_points_3.txt'/, 
     &               cfile4*18 /'ruler_points_4.txt'/ 

      offset1  = 12.36
      offset2  = -2.18
      offset3  = -7.036
      offset4  = 11.74
      scalesta = 10.42                                ! Matthew's scale for degree standard 
C       scalesta = 10.5466                              ! scale per degree standard from 0 to 56 degrees Bernie's

      open (13,file=cfile1,status='old',recl=120,access='sequential',form='formatted',iostat=iReadashi)
        if(iReadashi.ne.0) then
          print *, 'There is no ASHI ruler_points_1.txt file in this directory'
          close(13)
          stop
        end if
      NVAL = 0
      do I=1,INUM
        read (13,'(F13.0,F14.0)',iostat=iReadashi)  distx(1,I),disty(1,I)
        if(iReadashi.eq.0) then
          NVAL = NVAL + 1
        else
          close(13)
	  end if
      end do
      print *, ' '
 	print *, 'Number of ASHI ruler_points_1.txt read in =',NVAL

      open (13,file=cfile2,status='old',recl=120,access='sequential',form='formatted',iostat=iReadashi)
        if(iReadashi.ne.0) then
          print *, 'There is ASHI ruler_points_2.txt file in this directory'
          close(13)
          stop
        end if
      NVAL = 0
      do I=1,INUM
        read (13,'(F13.0,F14.0)',iostat=iReadashi)  distx(2,I),disty(2,I)
        if(iReadashi.eq.0) then
          NVAL = NVAL + 1
        else
          close(13)
	  end if    
      end do
 	print *, 'Number of ASHI ruler_points_2.txt read in =',NVAL

      open (13,file=cfile3,status='old',recl=120,access='sequential',form='formatted',iostat=iReadashi)
        if(iReadashi.ne.0) then
          print *, 'There is no ASHI ruler_points_3 file in this directory'
          close(13)
          stop
        end if
      NVAL = 0
      do I=1,INUM
        read (13,'(F13.0,F14.0)',iostat=iReadashi)  distx(3,I),disty(3,I)
        if(iReadashi.eq.0) then
          NVAL = NVAL + 1
        else
          close(13)
	  end if    
      end do
 	print *, 'Number of ASHI ruler_points_3.txt read in =',NVAL

      open (13,file=cfile4,status='old',recl=120,access='sequential',form='formatted',iostat=iReadashi)
        if(iReadashi.ne.0) then
          print *, 'There is no ASHI ruler_points_4.txt file in this directory'
          close(13)
          stop
        end if
      NVAL = 0
      do I=1,INUM
        read (13,'(F13.0,F14.0)',iostat=iReadashi)  distx(4,I),disty(4,I)
        if(iReadashi.eq.0) then
          NVAL = NVAL + 1
        else
          close(13)
	  end if    
      end do
 	print *, 'Number of ASHI ruler_points_4.txt read in =',NVAL

      do J=1,4
        do I=1,NVAL
          distx2 = (distx(J,I)-distx(J,1))**2
          disty2 = (disty(J,I)-disty(J,1))**2
          if(J.eq.1) then
            disttp1(I) = sqrt(distx2 + disty2) + offset1 
            degree1(I) = 0.900 * 5.0 * float(I-1) + offset1/scalesta
            if(I.eq.1) then
              print*,' '
              write(*,'(A,I1,A,F8.4,A,F7.3)')'#',J,'  DISTP     DiSTD   P/D ratio P/S*D ratio ratio =',scalesta,' offset1 =',offset1
            end if
            write(*,'(4F10.4)'),disttp1(I), degree1(I), disttp1(I)/degree1(I),(disttp1(I)/degree1(I))/scalesta
          end if
C          if(I.eq.NVAL) stop 
 
          if(J.eq.2) then
            disttp2(I) = sqrt(distx2 + disty2) + offset2
            degree2(I) = 0.900 * 5.0 * float(I-1) + offset2/scalesta             
            if(I.eq.1) then
              print*,' '
              write(*,'(A,I1,A,F8.4,A,F7.3)')'#',J,'  DISTP     DiSTD   P/D ratio P/S*D ratio ratio =',scalesta,' offset2 =',offset2
            end if
            write(*,'(4F10.4)'),disttp2(I),degree2(I),disttp2(I)/degree2(I), (disttp2(I)/degree2(I))/scalesta
          end if
          if(J.eq.3) then
            disttp3(I) = sqrt(distx2 + disty2) + offset3
            degree3(I) = 0.900 * 5.0 * float(I-1) + offset3/scalesta
            if(I.eq.1) then
              print*,' '
              write(*,'(A,I1,A,F8.4,A,F7.3)')'#',J,'  DISTP     DiSTD   P/D ratio P/S*D ratio ratio =',scalesta,' offset3 =',offset3
            end if
            write(*,'(4F10.4)'),disttp3(I),degree3(I),disttp3(I)/degree3(I), (disttp3(I)/degree3(I))/scalesta
          end if           
          if(J.eq.4) then
            disttp4(I) = sqrt(distx2 + disty2) + offset4 
            degree4(I) = 0.900 * 5.0 * float(I-1) + offset4/scalesta
            if(I.eq.1) then
              print*,' '
              write(*,'(A,I1,A,F8.4,A,F7.3)')'#',J,'  DISTP     DiSTD   P/D ratio P/S*D ratio ratio =',scalesta,' offset4 =',offset4
            end if
            write(*,'(4F10.4)'),disttp4(I),degree4(I),disttp4(I)/degree4(I), (disttp4(I)/degree4(I))/scalesta
          end if          
        end do
      end do

      do I=1,NVAL   
        if(I.eq.1) then
          print*,' '
          write(*,'(A,F8.4)')'Ave DISTP     DiSTD   P/D ratio P/S*D ratio  S*D/P    scale =',scalesta
        end if
        disttpa  = (disttp1(I) + disttp2(I) + disttp3(I) + disttp4(I))/4.0
        degreea  = (degree1(I) + degree2(I) + degree3(I) + degree4(I))/4.0      
        pixdegra = (disttp1(I)/degree1(I) + disttp2(I)/degree2(I) + disttp3(I)/degree3(I) + disttp4(I)/degree4(I))/4.0  
        pixdegpa = (disttp1(I)/degree1(I) + disttp2(I)/degree2(I) + disttp3(I)/degree3(I) + disttp4(I)/degree4(I))/(4.0*scalesta) 
        pdegpix  = 1.0/pixdegpa
        write(*,'(5F10.4)'),disttpa,degreea,pixdegra,pixdegpa,pdegpix
      end do
      stop
      end