*   -------------------------------------------------------------------
*     Written by  J.O'Brien  ( Department of Mathematics )  22/ 3/'88
*   -------------------------------------------------------------------
 
***********************************************************************
*    Sizes of the arrays specified.   These can be easily changed
*    to accomodate larger data files.

*       j1 :  size of arrays required to read all the data
*             from the input file

*       j2 :  size of arrays required to read data from the input
*             file corresponding to a particular set of coefficients

*       j3 :  size of arrays required to write estimates of
*             the coefficients to the output file  *.dat3

*       j4 :  corresponds to the maximum no. of N0 values used in
*             writing the coefficients to the output file  *.dat2

        parameter (j1 =2800) 
        parameter (j2 =150)
        parameter (j3 = 280)
        parameter (j4 = 159)

***********************************************************************

*    Those estimates that are acceptable for calculating the final 
*    estimates of the critical point and the critical exponent 
*    must have been calculated with more that a certain number of
*    values, jpn .

        parameter (jpn = 3)

*    Arrays set up and character variables defined

        real*8 cpdim(j1), cedim(j1),
     1         cpoint(j2), cexpon(j2), sr(j2), si(j2), er(j2), ei(j2),
     2         ranger, rangei, srkk, sikk, erkk, eikk,
     3         cp1, cp2, cp, eb, cal, cal1, cal2, factor, srr, check,
     4         xmean1(j3), ymean1(j3), sdx1(j3), sdy1(j3),
     5         xmean2(j3), ymean2(j3), sdx2(j3), sdy2(j3),
     6         sumx, sumx2, sumy, sumy2, exmin1, eymin1, exmin2, eymin2,
     7         qcr1, qcrer1, qex1, qexer1, qcr2, qcrer2, qex2, qexer2,
     8         pcr1, pcr2, pex1, pex2, ecr1, ecr2, eex1, eex2

        integer mount(3),fincr1(j3),fincr2(j3),finex1(j3),finex2(j3),
     1          na2(j1), na1(j1), na0(j1), nak(j1), nlabel(100),
     2          maxord(j1), lnum1(j3), lnum2(j3),
     3          ktotal,k1col1,k1col2,k1diff,k1tot

        character*78 title
        character*20 infile, outfi1, outfi2, outfi3
        character*15 out
        character*12 root2(j2), blank, double
        character*10 ex1, ey1, ex2, ey2, cpcon(j1), numdat(3,j4),
     1               xm1,xsd1,ym1,ysd1,xm2,xsd2,ym2,ysd2,
     2               mark(11), dbroot, flag04,
     3               pcr3, pcr4, pex3, pex4, ecr3, ecr4, eex3, eex4
        character*9 cecon(j1), mumdat(3,j4)
        character*8 flag01, flag02, loc01
        character*2 flag03, loc02
        character*1 ok,yes,no,answer,dr(3,j4),drdim(j1),drok,dd,star,
     1              lmark1(j3), lmark2(j3), recall, flag05, biased,
     2              first, ans1, ans2, ans3,flag15,flag16,ans15,ans16,
     3              ans99

*    Data initialized and flags defined

        data flag01, flag02, flag03 / '1 approx', '2 approx', '  '/
        data check /0.00000000001d+00/
        data n2, n1, n0, nk, iend, nodata /0,0,0,0,0,-99/
        data ok, yes, no, answer, recall / 'n', 'y', 'n', 'y', 'n'/

        data nlabel /1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
     1  21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,
     2  41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,
     3  61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,
     4  81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100/

        data mark /11*'----------'/
        data dbroot /'suspected '/
        data blank, double / '            ', ':double root'/
        data root2 /j2*'            '/
        data drok, star /' ', 'x'/
        data first /'y'/
        data flag15, flag16, ans15, ans16, ans99 /'n','n','y','y','n'/

        data eb,cal /0.0000001, 0.0000001/

        write (6,9876)
9876    format (/,1x,'NOTE: For accessing large input files, the size of'
     1  ,/,7x,'the arrays may have to be increased.')

*    Names of input and output files defined

        write (6,1001)
1001    format (//,' Name of input file : ',$)
        read (5,1002) infile
1002    format (a)
        write (6,1003)
1003    format (//,' Primary filename for result files : ',$)
        read (5,1002) out

        outfi1 = out//'.dat1'
        outfi2 = out//'.dat2'
        outfi3 = out//'.dat3'

        write (6,1074) 
1074    format (//,' Do you want the file  *.dat1  created (y/n) ? ',$)
        read (5,8001) ans3

*    Files opened

        open (31, file=infile, status='old', err=999)
        ok = yes
        if (ans3 .eq. yes) open (32, file=outfi1, status='unknown')
        open (33, file=outfi2, status='unknown')
        open (34, file=outfi3, status='unknown')

        read (31,1005) title
1005    format (1x,a)

*    Input file read until the first flag is located

        read (31,1006) loc01
1006    format (13x,a)
            do while ((loc01 .ne. flag01) .and. (loc01 .ne. flag02))
                read (31,1006) loc01
            end do

        backspace 31
        backspace 31

*    Check to see if it is a biased or unbiased case

        read (31,6501) flag05
6501    format (37x,a)
        if (flag05 .eq. 'D') then
            biased = yes
            backspace 31
            read (31,6502) cp
6502        format (4x,d37.0)
            write (6,6504) cp
6504        format (//,8x,6('>'),3x,'Biased Case :  Critical Point = ',
     1              f10.7,3x,6('<'),/)
        else
            biased = no
        end if

*    Values of the Critical Point and the search range defined

        if (biased .eq. no) then
            write (6,1004)
1004        format (//,' Input values of Critical Point and Range :',/)
            read (5,*) cp, eb
        end if

        write (6,8070)
8070    format (/,' Defective range factor for positive real axis : ',$)
        read (5,*) ranger
        write (6,8080)
8080    format (/,' Absolute defective range value for the complex plane
     1 : ',$)
        read (5,*) rangei
        write (6,1111)
1111    format (//,33x,'Please Wait',/)

            if (eb .lt. 0.0) then
                eb = 0.0
            end if

        cp1 = cp - eb
        cp2 = cp + eb

*    Title and Critical Point range read to output files

        if (ans3 .eq. yes) write (32,2001) title
        write (33,2002) title
2001    format (//,1x,a,/,1x,73('='),/)
2002    format (1x,a,/,1x,127('-'),/)

        if (biased .eq. no) then
          if (ans3 .eq. yes) write (32,2003) cp1, cp2
          write (33,2003) cp1, cp2
2003      format(1x,'Critical Point Range : ',f10.7,' <= Xc <= ',f10.7,/)
        else
          if (ans3 .eq. yes) write (32,6503) cp
          write (33,6503) cp
6503      format (1x,'Critical Point (biased) : ',f10.7,/)
        end if

        if (ans3 .eq. yes) write (32,3310) ranger, rangei
        write (33,3310) ranger, rangei
3310    format (1x,'Defective range factor for positive real axis : ',
     1     f7.3,/,1x,'Absolute defective range value for the complex
     2 plane : ',f11.7,/)
        if (ans3 .eq. yes) write (32,2005)
2005    format (1x,73('-'),/,3x,'N2',4x,'N1',4x,'N0',5x,'K',11x,
     1          'Critical Point',7x,'Critical Exponent',/,1x,73('-'),/)
        write (33,2006)
2006    format (1x,127('-'),/)

*     Check to see if data set is first or second order D.E.
        
        read (31,1007,iostat=iend) n0, n1, n2, nk
1007    format (47x,i2,2x,i2,2x,i2,2x,i2,//)
            if (loc01 .eq. flag01) then
                nk = n2
                n2 = 0
            end if

*    Read data until end of input file is encountered  

            loop = 0
            do while (iend .eq. 0)
                loop = loop + 1

                if (loc01 .eq. flag01) then
                    maxord(loop) = n0 + n1 + nk + 1
                else
                    maxord(loop) = n0 + n1 + n2 + nk + 2
                end if

*    Set counters equal to zero

                im = 0
                imax = 0

              do j = 1, n1
1414            read (31,1008) srkk, sikk, erkk, eikk, flag04
1008            format (4x,d19.0,1x,d19.0,1x,d18.0,1x,d18.0,4x,a)
                if (       (flag04 .eq. dbroot)
     1               .and. (abs(si(j-1)) .lt. check)    
     2               .and. (abs(ei(j-1)) .lt. check)
     3               .and. (sr(j-1) .ge. cp1)
     4               .and. (sr(j-1) .le. cp2)) then
                        root2(im) = double
                        read (31,8888) srr
8888                    format (4x,d19.0)
                        goto 1414
                else if (   flag04 .eq. dbroot) then
                        read (31,8888) srr
                        goto 1414
                else if (   (flag04 .ne. dbroot)
     1               .and. (abs(sikk) .lt. check)
     2               .and. (abs(eikk) .lt. check)
     3               .and. (srkk .ge. cp1)
     4               .and. (srkk .le. cp2)) then
                        sr(j) = srkk
                        si(j) = sikk
                        er(j) = erkk
                        ei(j) = eikk
                        im = im + 1
                        imax = im
                        cpoint(im) = sr(j)
                        cexpon(im) = er(j)
                        root2(im) = blank
                else 
                        sr(j) = srkk
                        si(j) = sikk
                        er(j) = erkk
                        ei(j) = eikk
                end if
              end do

*       Sort values of the Critical Point

                if (imax .gt. 1) then
                    call sortcp (cpoint, cexpon, root2, imax, cp)
                end if

                nak(loop) = nk
                na0(loop) = n0
                na1(loop) = n1
                na2(loop) = n2

*       Write the acceptable data to the output file

                if (imax .ne. 0) then

                call defect (sr,si,n1,cp,cp1,cp2,ranger,rangei,dd)

                drdim(loop) = dd

                    cpdim(loop) = cpoint(1)
                    cedim(loop) = cexpon(1)

                     nodata = 99

                if (ans3 .eq. yes) then
                   if (loc01 .eq. flag01) then
                   write(32,2007)n1,n0,nk,dd,cpoint(1),root2(1),cexpon(1)
2007            format(/,9x,i2,4x,i2,4x,i2,13x,a,1x,f10.7,1x,a,1x,f10.7)
                   else
                 write(32,2008)n2,n1,n0,nk,dd,cpoint(1),root2(1),cexpon(1)
2008        format(/,3x,i2,4x,i2,4x,i2,4x,i2,13x,a,1x,f10.7,1x,a,1x,f10.7)
                   end if
                end if

                if (ans3 .eq. yes) then
                   if (imax .gt. 1) then
                write(32,2009)(cpoint(i),root2(i),cexpon(i),i=2,imax)
2009                 format (38x,f10.7,1x,a,1x,f10.7)           
                   end if
                end if

                else

                    cpdim(loop) = 0.0
                    cedim(loop) = 0.0
                    drdim(loop) = drok

                if (ans3 .eq. yes) then
                   if (loc01 .eq. flag01) then
                       write (32,2010) n1, n0, nk
2010                   format (/,9x,i2,4x,i2,4x,i2,13x,'/',23x,'/')
                   else
                       write (32,2011) n2, n1, n0, nk
2011                   format (/,3x,i2,4x,i2,4x,i2,4x,i2,13x,'/',23x,'/')
                   end if
                end if

                end if

*       Read the next data set and check if first or second order D.E.

                read (31,1009,iostat=iend) loc01, n0, n1, n2, nk
1009            format (13x,a,26x,i2,2x,i2,2x,i2,2x,i2,//)

                    if (loc01 .eq. flag01) then
                        nk = n2
                        n2 = 0
                    end if

*    Return to primary loop   

            lmax = loop
            end do

        call conver (cpdim,cedim,lmax,cpcon,cecon)

*    Data read to file 33

*    Find the number of columns that are required for 1st order d.e.

        nn = 1
            do while (nak(nn) .eq. 1)
                nmax = na0(nn)
                nn = nn + 1
            end do
        ncol = nmax - na0(1) + 1

        if (ncol .gt. 11) then
            write (6,715)
715         format (/,1x,39('- '),/,
     1      22x,'Print Options For Output File *.Dat2',/,1x,39('- '),/)
        end if

        do while (ans15 .eq. yes)

704     if (ncol .gt. 11) then

            ans99 = yes
            write (6,701) na0(1), nmax, ncol
701         format (/,' First order D.E. has N0 from ',i2,' to ',i2,
     1      ' (',i2,')',/)
            write (6,702)
702         format(' Read in start and finish N0 values for printing :',/,
     1         ' Note :  total must be <= 11 at any one time',/)

            read (5,*) k1col1, k1col2

        else

            k1col1 = na0(1)
            k1col2 = nmax

        end if

        ktotal = k1col2 - k1col1 + 1
        k1diff = k1col1 - na0(1) + 1
        k1tot = k1col2 - na0(1) + 1

        if (ktotal .gt. 11) then
            write (6,703)
703     format(' The total no. of columns for printing must be <= 11')
            goto 704
        else if (k1col1 .ge. k1col2) then
            write (6,705)
705         format (' Error - input values again')
C**********Diagnostic
                write(6,7729)k1col1,k1col2
 7729           format(' k1col1=',i3,'  k1col2=',i3/)
C**************End diagnosticC********8
            goto 704
        else if ((k1col1 .lt. na0(1)) .or. (k1col2 .gt. nmax)) then
            write (6,705)
C**********Diagnostic
                write(6,7829)k1col1,k1col2,na0(1),nmax
 7829   format(' k1col1=',i3,'  k1col2=',i3,' na0(1)=',i3,' nmax=',i3/)
C**************End diagnosticC********8
            goto 704
        end if

        if (flag15 .eq. yes) write (33,2099)

        write (33,2012)
2012    format (1x,'FIRST  ORDER  DIFFERENTIAL  EQUATION  (K = 1)',/,
     1  1x,45('='),/)
        write (33,2013) (mark(i), i = 1, ktotal)
        write (33,2014) (nlabel(i), i = k1col1, k1col2)
        write (33,2013) (mark(i), i = 1, ktotal)
        write (33,2017)
2013    format (1x,'------',11(1x,a))
2014    format (1x,'N0 :',7x,11(i2,9x))

*    Next, find the highest value of k for the first order

        kk = 1
            do while ( (na2(kk) .eq. 0) .and. (kk .le. lmax) )
                kmax = nak(kk)
                kk = kk + 1
            end do

*    Next, print out data in tabular form for first order d.e.

            nend = 0
        do i = 1, kmax

        if (biased .eq. no) then
            if (i .eq. 5 .or. i .eq. 10 .or. i .eq. 15 .or. i .eq. 20)then
                write (33,2099)
2099            format ('1',//)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2014) (nlabel(k), k = k1col1, k1col2)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2017)
            end if
        else
            if(i .eq. 11 .or. i .eq. 21 .or. i .eq. 31 .or. i .eq. 41)then
                write (33,2099)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2014) (nlabel(k), k = k1col1, k1col2)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2017)
            end if
        end if

            write (33,2015) i
2015        format (1x,'L = ',i2)
            nstart = nend + 1
            n = nstart
            do while (nak(n) .eq. i)
                nend = n
                n = n + 1
            end do
                do j = 1, 3
                    do k = 1, ncol
                        nnn = j + (k - 1)*3 + nstart - 1
                            if (nnn .le. nend) then
                                mount (j) = k
                                numdat (j,k) = cpcon(nnn)
                                mumdat (j,k) = cecon(nnn)
                                dr(j,k) = drdim(nnn)
                            end if
                    end do
                end do
            do j = 1, 3
    
              if (mount(j) .gt. k1tot) mount(j) = k1tot

              if (biased .eq. no) then
                write (33,2016) (numdat(j,n), n=k1diff, mount(j))
                write (33,2216) (mumdat(j,n),dr(j,n), n=k1diff,mount(j))
                write (33,2017)
              else
                write (33,2216) (mumdat(j,n),dr(j,n), n=k1diff,mount(j))
              end if
2016            format (7x,11(1x,a))
2216            format (7x,11(1x,a,a))
2017            format (1x)
            end do
        end do

        if (ncol .gt. 11) then
            flag15 = yes
            write (6,711)
711     format(/,
     1  ' Do you wish to print more 1st Order results (y/n) : ',$)
            read (5,8001) ans15
        else
            ans15 = no
        end if

        end do

*    Number of data points for First Order

        lmax1 = nend

*    *******************************
        if ( lmax1 .ne. lmax ) then
*    *******************************

*    Find the number of columns that are required for 2nd order

        nend1 = nend + 1
        nn = nend + 1
            do while (nak(nn) .eq. 1)
                nmax = na0(nn)
                nn = nn + 1
            end do
        ncol = nmax - na0(nend1) + 1

        do while (ans16 .eq. yes)

707     if (ncol .gt. 11) then

            ans99 = yes
            write (6,706) na0(nend1), nmax, ncol
706         format (/,' Second order D.E. has N0 from ',i2,' to ',i2,
     1      ' (',i2,')',/)

            write (6,702)

            read (5,*) k1col1, k1col2

        else

            k1col1 = na0(nend1)
            k1col2 = nmax

        end if

        ktotal = k1col2 - k1col1 + 1
        k1diff = k1col1 - na0(nend1) + 1
        k1tot = k1col2 - na0(nend1) + 1

        if (ktotal .gt. 11) then
            write (6,703)
            goto 707
        else if (k1col1 .ge. k1col2) then
            write (6,705)
            goto 707
        else if ((k1col1 .lt. na0(nend1)) .or. (k1col2 .gt. nmax)) then
            write (6,705)
            goto 707
        end if

        write (33,2018)
2018    format ('1',//,1x,'SECOND  ORDER  DIFFERENTIAL  EQUATION  (K = 2)'
     1  ,/,1x,46('='),/)
        write (33,2013) (mark(i), i = 1, ktotal)
        write (33,2014) (nlabel(i), i = k1col1, k1col2)
        write (33,2013) (mark(i), i = 1, ktotal)
        write (33,2017)

*    Next, find the highest value of k for the second order d.e.

                kmax = nak(lmax)

*    Next, print out data in tabular form for second order d.e.

            nend = nend1 - 1
        do i = 1, kmax

        if (biased .eq. no) then
            if (i .eq. 5 .or. i .eq. 10 .or. i .eq. 15 .or. i .eq. 20)then
                write (33,2099)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2014) (nlabel(k), k = k1col1, k1col2)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2017)
            end if
        else
            if(i .eq. 11 .or. i .eq. 21 .or. i .eq. 31 .or. i .eq. 41)then
                write (33,2099)
                write (33,2013) (mark(k), k = 1, ktotal )
                write (33,2014) (nlabel(k), k = k1col1, k1col2)
                write (33,2013) (mark(k), k = 1, ktotal)
                write (33,2017)
            end if
        end if

            write (33,2015) i
            nstart = nend + 1
            n = nstart
            do while ((nak(n) .eq. i) .and. (n .le. lmax))
                nend = n
                n = n + 1
            end do
                do j = 1, 3
                    do k = 1, ncol
                         nnn = j + (k - 1)*3 + nstart - 1
                            if (nnn .le. nend) then
                                mount (j) = k
                                numdat (j,k) = cpcon(nnn)
                                mumdat (j,k) = cecon(nnn)
                                dr (j,k) = drdim(nnn)
                            end if
                    end do
                end do

            if ((na1(nend).eq.na2(nend-1)).and.(nak(nend).eq.nak(nend-1)))
     1        then
                mount (1) = mount (1) + 1
                numdat (1,mount(1)) = numdat (3,mount (3))
                mumdat (1,mount(1)) = mumdat (3,mount (3))
                mount (3) = mount (3) - 1
            else if (nak(nend) .ne. nak(nend - 1)) then
                mount (2) = 0
                mount (3) = 0
            end if

            do j = 1, 3

            if (mount(j) .gt. k1tot) mount(j) = k1tot

              if (biased .eq. no) then
                write (33,2016) (numdat(j,n),n=k1diff,mount(j))
                write (33,2216) (mumdat(j,n),dr(j,n),n=k1diff,mount(j))
                write (33,2017)
              else
                write (33,2216) (mumdat(j,n),dr(j,n),n=k1diff,mount(j))
              end if
            end do
        end do

        if (ncol .gt. 11) then
            flag16 = yes
            write (6,713)
713     format(/,
     1  ' Do you wish to print more 2nd Order results (y/n) : ',$)
            read (5,8001) ans16
        else
            ans16 = no
        end if

        end do

        if (ans99 .eq. yes) then
            write (6,716)
716         format (/,1x,39('- '),/,
     1      19x,'End Of Print Options For Output File *.Dat2',/,
     2      1x,39('- '),/)
        end if

*    ************************
        end if
*    ************************

*    Calculation of estimates of Critical Point and Exponent

        jormin = 9999
        jormax = -9999
        do i = 1, lmax
            if (maxord(i) .lt. jormin) jormin = maxord(i)
            if (maxord(i) .gt. jormax) jormax = maxord(i)
        end do

        if (biased .eq. no) then
            write (6,2999)
2999        format (/,1x,'Error range for accepting Xc estimates : ')
            read (5,*) cal
        end if

        do while (answer .eq. yes)

        if ((recall .eq. yes) .and. (biased .eq. no)) then
            write (6,3000) cal
3000        format (/,1x,'Error range for accepting Xc estimates : ',/,
     1      1x,'( Range for previous acceptance was ',f10.7,' )',/)
            read (5,*) cal
        end if
        recall = yes

        cal1 = cp - cal
        cal2 = cp + cal

*    Calculate the means and twice the standard deviations for the
*    Critical Point and the Critical Exponent for different maximal
*    orders of the coefficients used in the construction of the
*    approximants.

        lmin2 = lmax1 + 1
        do i = jormin, jormax
            lnum = 0
            sumx = 0.0
            sumx2 = 0.0
            sumy = 0.0
            sumy2 = 0.0
                do j = 1, lmax1
                    if ((maxord(j).eq.i).and.(cpdim(j).ge.cal1)
     1                .and.(cpdim(j).le.cal2).and.(drdim(j).eq.drok))then
                        lnum = lnum + 1
                        sumx = sumx + cpdim(j)
                        sumx2 = sumx2 + cpdim(j)**2
                        sumy = sumy + cedim(j)
                        sumy2 = sumy2 + cedim(j)**2
                    end if
                end do
            if (lnum .ne. 0) then
              xmean1(i) = sumx/lnum
              ymean1(i) = sumy/lnum
              if (lnum .gt. 1) then
                sdx1(i) = (sumx2 - lnum*xmean1(i)**2)/(lnum - 1)
                sdx1(i) = 2.0*(abs(sdx1(i)))**0.5
                sdy1(i) = (sumy2 - lnum*ymean1(i)**2)/(lnum - 1)
                sdy1(i) = 2.0*(abs(sdy1(i)))**0.5
              else
                sdx1(i) = 0.0
                sdy1(i) = 0.0
              end if
            else 
              xmean1(i) = 0.0
              ymean1(i) = 0.0
              sdx1(i) = 0.0
              sdy1(i) = 0.0
            end if
            lnum1(i) = lnum

            lnum = 0
            sumx = 0.0
            sumx2 = 0.0
            sumy = 0.0
            sumy2 = 0.0

                do j = lmin2, lmax
                    if ((maxord(j).eq.i).and.(cpdim(j).ge.cal1)
     1                .and.(cpdim(j).le.cal2).and.(drdim(j).eq.drok))then
                        lnum = lnum + 1
                        sumx = sumx + cpdim(j)
                        sumx2 = sumx2 + cpdim(j)**2
                        sumy = sumy + cedim(j)
                        sumy2 = sumy2 + cedim(j)**2
                    end if
                end do
            if (lnum .ne. 0) then
              xmean2(i) = sumx/lnum
              ymean2(i) = sumy/lnum
              if (lnum .gt. 1) then
                sdx2(i) = (sumx2 - lnum*xmean2(i)**2)/(lnum - 1)
                sdx2(i) = 2.0*(abs(sdx2(i)))**0.5
                sdy2(i) = (sumy2 - lnum*ymean2(i)**2)/(lnum - 1)
                sdy2(i) = 2.0*(abs(sdy2(i)))**0.5
              else
                sdx2(i) = 0.0
                sdy2(i) = 0.0
              end if
            else
              xmean2(i) = 0.0
              ymean2(i) = 0.0
              sdx2(i) = 0.0
              sdy2(i) = 0.0
            end if
            lnum2(i) = lnum
        end do

        do i = jormin, jormax
            if (lnum1(i) .le. jpn) then
                lmark1(i) = star
            else
                lmark1(i) = drok
            end if
            if (lnum2(i) .le. jpn) then
                lmark2(i) = star
            else
                lmark2(i) = drok
            end if
        end do

*    Calculate final estimates of the Critical Point and the Exponent

        exmin1 = 9999.
        eymin1 = 9999.
        exmin2 = 9999.
        eymin2 = 9999.

        ans1 = no
        ans2 = no
        do i = jormin, jormax
          if((sdx1(i).lt.exmin1).and.(lnum1(i).gt.jpn)) then
              exmin1 = sdx1(i)
              ans1 = yes
          end if
          if((sdy1(i).lt.eymin1).and.(lnum1(i).gt.jpn)) eymin1=sdy1(i)
          if((sdx2(i).lt.exmin2).and.(lnum2(i).gt.jpn)) then
              exmin2 = sdx2(i)
              ans2 = yes
          end if
          if((sdy2(i).lt.eymin2).and.(lnum2(i).gt.jpn)) eymin2=sdy2(i)
        end do

        if ((ans1.eq.no).and.(ans2.eq.no)) then
            goto 7777
        end if

        if (first .eq. yes) then
        write (6,5001)
5001    format (/,1x,'The minimum errors are as follows : ',/)

        if (biased .eq. no) then
        if (ans1 .ne. no) then
        write (6,5002) exmin1, eymin1
5002    format(/,1x,'K = 1',10x,'Minimum Critical Point Error : ',f10.7,
     1         /,1x,'K = 1',10x,'Minimum Exponent Error       : ',f10.7)
        end if
        if (ans2 .ne. no) then
        write (6,5052) exmin2, eymin2
5052    format(/,1x,'K = 2',10x,'Minimum Critical Point Error : ',f10.7,
     1         /,1x,'K = 2',10x,'Minimum Exponent Error       : ',f10.7)
        end if

        else

        if (ans1 .ne. no) then
        write (6,6505) eymin1
6505    format(/,1x,'K = 1',10x,'Minimum Exponent Error : ',f10.7,/)
        end if
        if (ans2 .ne. no) then
        write (6,6535) eymin2
6535    format(/,1x,'K = 2',10x,'Minimum Exponent Error : ',f10.7,/)
        end if

        end if

        end if
        first = no

        write (6,5003)
5003    format (/,1x,'Read in the error acceptance factor : ',
     1    /,1x,'( Resulting acceptance range = min error x factor )',/)
        read (5,5004)factor
5004    format (f10.0)

        exmin1 = exmin1*factor
        eymin1 = eymin1*factor
        exmin2 = exmin2*factor
        eymin2 = eymin2*factor

        qcr1 = 0.0
        qcrer1 = 0.0
        qex1 = 0.0
        qexer1 = 0.0
        qcr2 = 0.0
        qcrer2 = 0.0
        qex2 = 0.0
        qexer2 = 0.0
        mx1 = 0
        my1 = 0
        mx2 = 0
        my2 = 0

        do i = jormin, jormax

            if ((sdx1(i) .le. exmin1).and.(lnum1(i) .gt. jpn)
     1          .and. (sdx1(i) .gt. check)) then
                mx1 = mx1 + 1
                qcr1 = qcr1 + xmean1(i)/sdx1(i)
                qcrer1 = qcrer1 + 1.0/sdx1(i)
                fincr1(mx1) = i
            end if

            if ((sdx2(i) .le. exmin2).and.(lnum2(i) .gt. jpn)
     1          .and. (sdx2(i) .gt. check)) then
                mx2 = mx2 + 1
                qcr2 = qcr2 + xmean2(i)/sdx2(i)
                qcrer2 = qcrer2 + 1.0/sdx2(i)
                fincr2(mx2) = i
            end if

            if ((sdy1(i) .le. eymin1).and.(lnum1(i) .gt. jpn)
     1          .and. (sdy1(i) .gt. check)) then
                my1 = my1 + 1
                qex1 = qex1 + ymean1(i)/sdy1(i)
                qexer1 = qexer1 + 1.0/sdy1(i)
                finex1(my1) = i
            end if

            if ((sdy2(i) .le. eymin2).and.(lnum2(i) .gt. jpn)
     1          .and. (sdy2(i) .gt. check)) then
                my2 = my2 + 1
                qex2 = qex2 + ymean2(i)/sdy2(i)
                qexer2 = qexer2 + 1.0/sdy2(i)
                finex2(my2) = i
            end if
            
        end do

        if (mx1 .ne. 0) then
            pcr1 = qcr1/qcrer1
            ecr1 = (mx1**0.5)/qcrer1
        end if

        if (mx2 .ne. 0) then
            pcr2 = qcr2/qcrer2
            ecr2 = (mx2**0.5)/qcrer2
        end if

        if (my1 .ne. 0) then
            pex1 = qex1/qexer1
            eex1 = (my1**0.5)/qexer1
        end if

        if (my2 .ne. 0) then
            pex2 = qex2/qexer2
            eex2 = (my2**0.5)/qexer2
        end if

*    Results sent to file 34

7777    continue

        exmin1 = exmin1/factor
        eymin1 = eymin1/factor
        exmin2 = exmin2/factor
        eymin2 = eymin2/factor

        call chang1 (exmin1,eymin1,exmin2,eymin2,ans1,ans2,
     1               ex1,ey1,ex2,ey2)

        call chang2 (pcr1,ecr1,pex1,eex1,pcr2,ecr2,pex2,eex2,ans1,ans2,
     1               pcr3,ecr3,pex3,eex3,pcr4,ecr4,pex4,eex4)

        if (biased .eq. no) then

3333    format (1x,'|',4x,'|',56x,'|',4x,'|',56x,'|')
4444    format (1x,'|',61x,'|',61x,'|')
5555    format (1x,125('-'))
        write (34,3001) title
3001    format ('1',//,1x,a,/,1x,125('-'),/)
        write (34,3002) cal1, cal2
3002    format (1x,'Critical Point range used for calculations : ',
     1  f10.7,' <= Xc <= ',f10.7)
        write (34,4800) cp1, cp2
4800    format (1x,'(range used for initial data extraction    : ',
     1  f10.7,' <= Xc <= ',f10.7,')',/)
        write (34,3310) ranger, rangei
        write (34,3320)
3320    format (1x,125('-'))
        write (34,3003)
3003    format (1x,'|',28x,'K = 1',28x,'|',28x,'K = 2',28x,'|')
        write (34,5555)
        write (34,3004)
3004    format (1x,'|',2x,'n',1x,'|',6x,'Critical Point',10x,
     1  'Critical Exponent',6x,'L',2x,'|',2x,'n',1x,'|',6x,
     2  'Critical Point',10x,'Critical Exponent',6x,'L',2x,'|')
        write (34,5555)
        write (34,3333)

        do i = jormin, jormax

        call chang3 (xmean1(i),sdx1(i),ymean1(i),sdy1(i),lnum1(i),
     1       xmean2(i),sdx2(i),ymean2(i),sdy2(i),lnum2(i),
     2       xm1,xsd1,ym1,ysd1,xm2,xsd2,ym2,ysd2)

          write(34,3005)i,xm1,xsd1,ym1,ysd1,lnum1(i),lmark1(i),
     1                  i,xm2,xsd2,ym2,ysd2,lnum2(i),lmark2(i)
3005      format(1x,'|',1x,i2,1x,'|',2x,a,1x,a,5x,a,1x,a,3x,i2,a,1x,
     1              '|',1x,i2,1x,'|',2x,a,1x,a,5x,a,1x,a,3x,i2,a,1x,'|')

        end do

        write (34,3333)
        write (34,5555)

        if ((ans1 .eq. no) .and. (ans2 .eq. no)) then
            goto 1122
        end if

        write (34,7001) ex1, ey1, ex2, ey2
7001    format (3x,'Minimum Error :',2x,a,16x,a,26x,a,16x,a)
        write (34,3007)
3007    format (1x,125('-'),/,56x,'Final Estimates')
        write (34,3017) factor
3017    format (47x,'(error acceptance factor : ',f4.1,')',/)

        write (34,3008) pcr3, ecr3, pex3, eex3, pcr4, ecr4, pex4, eex4
3008    format (9x,a,1x,a,5x,a,1x,a,15x,a,1x,a,5x,a,1x,a)

3009    format (1x,125('-'),/)
        write (34,4600)
4600    format (/,1x,125('-'))

        write (34,4500)
4500    format(1x,'Values(n) used in calculating the final estimates :',/)

        if (mx1 .ne. 0) then
            write (34,4100) (fincr1(i), i = 1, mx1)
        end if

        if (my1 .ne. 0) then
            write (34,4200) (finex1(i), i = 1, my1)
        end if

        write (34,2017)

        if (mx2 .ne. 0) then
            write (34,4300) (fincr2(i), i = 1, mx2)
        end if

        if (my2 .ne. 0) then
            write (34,4400) (finex2(i), i = 1, my2)
        end if

4100    format (1x,'Critical Point    ( K = 1 ) : ',30(1x,i2,','))
4200    format (1x,'Critical Exponent ( K = 1 ) : ',30(1x,i2,','))
4300    format (1x,'Critical Point    ( K = 2 ) : ',30(1x,i2,','))
4400    format   (1x,'Critical Exponent ( K = 2 ) : ',30(1x,i2,','))
        write (34,3009)

        else

        write (34,6506) title
6506    format ('1',//,1x,a,/,1x,79('-'),/)
        write (34,6507) cp
6507    format (1x,'Critical Point (biased) : ',f10.7,/)
        write (34,3310) ranger, rangei
        write (34,6508)
6508    format (1x,79('-'))
        write (34,6509)
6509    format (1x,'|',17x,'K = 1',16x,'|',17x,'K = 2',16x,'|')
        write (34,6508)
        write (34,6510)
6510    format (1x,'|',2x,'n',1x,'|',6x,'Critical Exponent',6x,'L',3x,
     1  '|',2x,'n',1x,'|',6x,'Critical Exponent',6x,'L',3x,'|')
        write (34,6508)
        write (34,6511)
6511    format (1x,'|',4x,'|',33x,'|',4x,'|',33x,'|')
        do i = jormin, jormax

        call chang4 (ymean1(i),sdy1(i),lnum1(i),ymean2(i),sdy2(i),
     1               lnum2(i),ym1,ysd1,ym2,ysd2)

            write (34,6512) i,ym1,ysd1,lnum1(i),lmark1(i),
     1                      i,ym2,ysd2,lnum2(i),lmark2(i)
6512        format (1x,'|',1x,i2,1x,'|',3x,a,2x,a,3x,i2,a,
     1              2x,'|',1x,i2,1x,'|',3x,a,2x,a,3x,i2,a,2x,'|')
        end do
        write (34,6511)
        write (34,6508)

        if ((ans1 .eq. no) .and. (ans2 .eq. no)) then
            goto 1122
        end if

        write (34,6513) ey1, ey2
6513    format (3x,'Minimum Error :',4x,a,29x,a)
        write (34,6508)
        write (34,6514) factor
6514    format (33x,'Final Estimates',/,
     1          24x,'(error acceptance factor : ',f4.1,')',/)

        write (34,6515) pex3, eex3, pex4, eex4
6515    format (10x,a,2x,a,17x,a,2x,a,/)

        write (34,6508)
        write (34,6516)
6516    format
     1  (1x,'Values(n) used in calculating the final estimates :'/)

        if (my1 .ne. 0) then
            write (34,4200) (finex1(i), i = 1, my1)
        end if

        write (34,2017)

        if (my2 .ne. 0) then
            write (34,4400) (finex2(i), i = 1, my2)
        end if

        write (34,6508)

        end if

*    Perform calculations with a different range ?

        write (6,8000)
8000    format (/,1x,'Do you wish to re-estimate the parameters (y/n) ? '
     1  ,$)
        read (5,8001) answer
8001    format(a)

        end do

*    Go To statement number 1122

1122    continue

*    Check if any data was accepted

        if (nodata .lt. 0) then
            write (6,2019)
            if (ans3 .eq. yes) write (32,2019)
            write (33,2019)
2019        format (/,1x,73('-'),//,
     1              20x,'N O     D A T A     A C C E P T E D',/)
        else
            write (6,2020)
            if (ans3 .eq. yes) write (32,2020)
            write (33,2021)
2020        format (/,1x,73('-'),//,25x,'E N D     O F     D A T A',/)
2021        format (/,1x,127('-'),//,50x,'E N D     O F     D A T A',/)
        end if

*    Check if an error occurred in accessing input file

999      if (ok .eq. no) then
            write (6,9999)      
9999        format (//' Warning: Error occurred in opening input file -'
     1               /'           is the filename correct ?'//
     2                 ' Execution terminated'/)
            close (31)
            if (ans3 .eq. yes) close (32)
            close (33)
            close (34)
            stop
        end if

*    Close input data files

        close (31)
        if (ans3 .eq. yes) close (32)
        close (33)
        close (34)

*    Program terminated

        stop '                         E N D     O F     R U N'
        end

*    Subroutine to sort values of Critical Point

        subroutine sortcp (value1, value2, root2, numval, cp)

        real*8 value1 (numval), value2 (numval), diff (60)
        real*8 cp, temp, temp1, temp2
        character*12 root2 (numval), temp3
        integer flag

        do i = 1, numval
                diff(i) = abs(value1(i) - cp)
        end do
                
        listln = numval
        do while (listln .gt. 1)
            listln = listln/2
            flag = 1
                do while (flag .ne. 0)
                    flag = 0
                    numcom = numval - listln
                        do j=1,numcom
                            if (diff(j) .gt. diff(j + listln)) then
                                temp = diff(j)
                                temp1 = value1(j)
                                temp2 = value2(j)
                                temp3 = root2(j)
                                diff(j) = diff(j + listln)
                                value1(j) = value1(j + listln)
                                value2(j) = value2(j + listln)
                                root2(j) = root2(j + listln)
                                diff(j + listln) = temp
                                value1(j + listln) = temp1
                                value2(j + listln) = temp2
                                root2(j + listln) = temp3
                                flag = 1
                            end if
                        end do
                end do
        end do
        return
        end


*    This subroutine determines whether a value of the Critical Point
*    is defective.

        subroutine defect(sr,si,n1,cp,cp1,cp2,ranger,rangei,flag)

        real*8 sr(n1), si(n1), cp, cp1, cp2, ranger, rangei, upper
        character*1 flag, ok, notok

        data ok, notok /' ', '*'/

        flag = ok

        upper = ranger*cp

        do i = 1, n1

            if((((sr(i).lt.cp1).and.(sr(i).gt.0.0)).or.((sr(i).gt.cp2)
     1         .and.(sr(i).le.upper))) .and. (abs(si(i)).le.rangei))
     2      then

                flag = notok
            end if
        end do
        return
        end
    
*    Program to convert from numeric to character data
    
        subroutine conver (cpdim,cedim,lmax,cpcon,cecon)

        real*8 cpdim (lmax), cedim (lmax)
        character*10 cpcon (lmax)
        character*9 cecon (lmax)
        do i = 1, lmax
            write (cpcon(i),1000) cpdim(i)
            write (cecon(i),1001) cedim(i)
1000        format (f10.7)
1001        format (f9.6)
        end do
        do i = 1, lmax
          if (cpcon(i).eq.' 0.0000000'.and.cecon(i).eq.' 0.000000')then
              cpcon(i) = '   ----   '
              cecon(i) = '   ----  '
          end if
        end do
        return
        end


        subroutine chang1 (exmin1,eymin1,exmin2,eymin2,ans1,ans2,
     1  ex1,ey1,ex2,ey2)

        real*8 exmin1, eymin1, exmin2, eymin2
        character*10 ex1, ey1, ex2, ey2, none
        character*1 ans1, ans2, no
        data none /'     -    '/
        data no /'n'/

        if (ans1 .ne. no) then
            write (ex1,1000) exmin1
            write (ey1,1000) eymin1
        else
            ex1 = none
            ey1 = none
        end if
        if (ans2 .ne. no) then
            write (ex2,1000) exmin2
            write (ey2,1000) eymin2
        else
            ex2 = none
            ey2 = none
        end if
1000    format (f10.7)

        return
        end


        subroutine chang2 (pcr1,ecr1,pex1,eex1,pcr2,ecr2,pex2,eex2,
     1  ans1,ans2,pcr3,ecr3,pex3,eex3,pcr4,ecr4,pex4,eex4)

        real*8       pcr1,ecr1,pex1,eex1,pcr2,ecr2,pex2,eex2
        character*10 pcr3,ecr3,pex3,eex3,pcr4,ecr4,pex4,eex4,none
        character*1  ans1, ans2, no
        data none /'     -    '/
        data no /'n'/

        if (ans1 .ne. no) then
            write (pcr3,1000) pcr1
            write (ecr3,1000) ecr1
            write (pex3,1000) pex1
            write (eex3,1000) eex1
        else
            pcr3 = none
            ecr3 = none
            pex3 = none
            eex3 = none
        end if

        if (ans2 .ne. no) then
            write (pcr4,1000) pcr2
            write (ecr4,1000) ecr2
            write (pex4,1000) pex2
            write (eex4,1000) eex2
        else
            pcr4 = none
            ecr4 = none
            pex4 = none
            eex4 = none
        end if
1000    format (f10.7)

        return
        end


        subroutine chang3 (xxm1,xxsd1,yym1,yysd1,lnum1,
     1                     xxm2,xxsd2,yym2,yysd2,lnum2,
     2                     xm1,xsd1,ym1,ysd1,xm2,xsd2,ym2,ysd2)

        real*8 xxm1,xxsd1,yym1,yysd1,xxm2,xxsd2,yym2,yysd2
        integer lnum1, lnum2
        character*10 xm1,xsd1,ym1,ysd1,xm2,xsd2,ym2,ysd2,blank
        data blank /'     -    '/

        write (xm1,1000) xxm1
        write (xsd1,1000) xxsd1
        write (ym1,1000) yym1
        write (ysd1,1000) yysd1
        write (xm2,1000) xxm2
        write (xsd2,1000) xxsd2
        write (ym2,1000) yym2
        write (ysd2,1000) yysd2
1000    format (f10.7)

        if (lnum1 .eq. 0) then
            xm1 = blank
            xsd1 = blank
            ym1 = blank
            ysd1 = blank
        else if (lnum1 .eq. 1) then
            xsd1 = blank
            ysd1 = blank
        end if

        if (lnum2 .eq. 0) then
            xm2 = blank
            xsd2 = blank
            ym2 = blank
            ysd2 = blank
        else if (lnum2 .eq. 1) then
            xsd2 = blank
            ysd2 = blank
        end if

        return
        end


        subroutine chang4 (yym1,yysd1,lnum1,yym2,yysd2,lnum2,
     1                     ym1,ysd1,ym2,ysd2)

        real*8 yym1,yysd1,yym2,yysd2
        integer lnum1, lnum2
        character*10 ym1,ysd1,ym2,ysd2,blank
        data blank /'     -    '/

        write (ym1,1000) yym1
        write (ysd1,1000) yysd1
        write (ym2,1000) yym2
        write (ysd2,1000) yysd2
1000    format (f10.7)

        if (lnum1 .eq. 0) then
            ym1 = blank
            ysd1 = blank
        else if (lnum1 .eq. 1) then
            ysd1 = blank
        end if

        if (lnum2 .eq. 0) then
            ym2 = blank
            ysd2 = blank
        else if (lnum2 .eq. 1) then
            ysd2 = blank
        end if

        return
        end