! fatigue crack growth rates programme
! this programme file name--GB. FUR
! m——total test points
! Pmax——maximun load
! Pmin——minimun load
! b—— specimen thickness
! w——specimen widthness
! a0——initial crack length
! ys——yield strength
! Sratio—— Stress ratio
! Ni(l)——load cycle number
! Ai(1)—— a set of crack length
        dimension ni(100),ai(100),dadn(100)
        dimension delk(100),aa(10),nn(10),bb(3),name(20)
        integer qq;
        parameter (PI=3.1415926)
        Character(15) finnante,foutname,key5,ancurve15,dkdadn15,name50
        write(*,*)'input Title:'
        read(*, '(A)')name(8)
        write(*,*)'input data file name:'
        read(*,'(A)')finname
        write(*,*)'output data file name'
        read(*,'(A)')foutname
        write( *,*)'enter A-N curve name;'
        read(*,'(A)')ancurye
        write(*,*)'enter DELT K-da/dN curve name:'
        read(*,'(A)')dkdadn
       open (15,file=finname,status='old')
       open (16,file=foutname,status='new')
       open (17,file=ancurve,status='new')
       open (18,file=dkdadn,status='new')
       write(*,*)'Material Name;'
       read( *, '(A)')name(1)
       write(*,*)'No.of Specimen'
       read(*,'(A)')name(2)
       write(*,*)'environment'
       read(*,'(A)')name(3)
       write(*,*)'test frequency'
       read(*,'(A)')name(9)
       write(*,*)'Temperature'
       read(*,'(A)')name(4)
       write(*,*)'Moisture'
       read (*,'(A)')name(5)
       write(*,*)'Test Date'
       read(*,'(A)')name(6)
       write (*,*)'Testers'
       read(*,'(A)')name(7)
       write(16,300)name (8)
       format(lx,'******TEST REPORT:',A35,'******')
       write(16,299)
       format('/')
       write (16,400)name(1)
       format(1X,'Material Name&Orientation:',2x,A40)
       write(16,410) name(2)
       format(lx,'No.of Specimcn:',2x,A20)
       write(16,420) name(3)
       format (1x,'environments:',2x,A20)
       write(16,430)name(9)
       format(lx,'Test Frequency:',3x,A3,'Hz')
       write(l6,440)name(4)
       format(lx,'Temperature:',2x,A3,'c' )
       write(16,450)name(5)
       format (lx,'Moisture:',2x,A4,'%','RH')
       write(16,460)name(6)
       format(lx,'Test Date:',2x,A15)
       write(16,470)name (7)
       format (lx,'Testers:',2x,A20)
       write(*,*)'asking for specimen kind:SEB,CT OR MT'
       read(*,'(A)')key
! begin to operate seven points poly
       read(15,*)m,pmax.pmin,b,w,ys,a0,Sratio
       write(16,20)
       format(/)
       write(16,21)
       format(5x,'=========',lx,'RESULTS OF TEST DATA PROCESSING',&
       1x,'========')
       write(16,22)
       format(/)
       write(16,25)m,pmax,pmin,b,w,a0, ys, Stratio
       format(lx,'the test points:',I4,2x,'pmax=',f8.1,'N', 2x,'pmin=',&
       f8.1.'N',/,1x,'b=',f5.2,'mm',2x,'w=', f6.2,'mm',2x,'A0=',f6.3,'mm',&
       /lx,'yield limit ys=',f8.2,'MPa',5x,'Stress Ratio R=',f6.3)
       read(15,*)(Ni(I),Ai(I),I=1,m)
       DO 30 I=1,m
       AI(I)=A0+AI(I)
       Continue
       write(16,105)
       format(lx,'cbs,no',lx,'cycle',2x,'a(means)',2x,'a(areg)',3x,'m.c.c.',4x,'delk',6x,'dadn')
       k=o
       pp=pmax-pmm
       r=pmin/pmax
       write(16,95)(i,Ni(i),Ai(1),I=1,3)
       m=m-6
       do 100 i=1,m
       i=0
       k=k+1
      kl=k+6
       do 64 j=k,k1
       i=i+1
       aa(1)=ai(j)
       nn(1)=ni(j)
       continue
       c1=0.5*(nn(1)+nn(7))
       c2=0.5*(nn(7)-nn(1))
       sx= 0.0
       sx2=0.0
       sx3=0.0
       sx4=0.0
       sy=0.0
       syx=0.0
       syx2=0.0
       do 70 j=1,7
       x=(NN(j)-cl)/c2
       yy=aa(j)
       sx=sx+x
       sx2=sx2+x**2
       sx3=sx3+x**3
       sx4=sx4+x**4
       sy=sy+yy
       syx=syx+x*yy
       syx2=syx2+yy*(x**2)
       continue
       den=7.0*(sx2*sx4-sx3**2)-sx*(sx*sx4-sx2*sx3)+sx2*(sx*sx3-sx2**2)
       t2=sy*(sx2*sx4-sx3**2)-syx*(sx*sx4-sx2*sx3)+syx2*(sx*sx3-sx2**2)
       bb(1)=t2/den
       t3=7.0*(syx*sx4-syx2*sx3)-sx*(sy*sx4-syx2*sx2)+sx2*(sy*sx3-syx*sx2)
       bb(2)=t3/den
       t4=7.0*(sx2*syx2-sx3*syx)-sx*(sx*syx2-sx3*sy)+sx2*(sx*syx-sx2*sy)
       bb(3)=t4/den
       yb=sy/7.0
       rss=0.0
       tss=0.0
       do 75 j=1,7
       x=(nn(j)-cl)/c2
       yhat=bb(1)+bb(2)*x+bb(3)*(x**2)
       rss=rss+(aa(j)-yhat)**2
       tss =tss+(aa(j)-yb)**2
       continue
       r2=1.0-rss/tSS
       dadn(i)=(bb(1)/c2+2.0*bb(2)*(nn(4)-c1)/(c2**2))
       x=(nn(4)-cl)/c2
       ar=bb(1)+bb(2)*x+bb(3)*(x**2)
       s=1e10
       snet=0.0
       qq=i+3
       IF(key.se.'SEB'. or. key. eq.'seb' )THEN
       t=ar/w
       ft=6*t**0.5/(1.0+2.0*t)/(1.0-t)**1.5*(1.99-t*(1-t)*(2.15-3.93*t)/(2.7*t**2))
       delk(i)=(ft*pp((b*sqrt(w))))/(10**1.5)
       s=ys*(w*(1-t))** 2
       ax=pmax*6.0*w/b
       snet=pmax/(b*w*(l-t))
       else
       if(key.eq.'ct'.or.key.eq.'CT')then
       t=ar/w
       ft=(2+t)*(0.8866+4.64*t-13.32*t*t+14.72*t**3-5.6*t**4)/((1-t)**1.5)
       delk(i)=(ft*pp)/(b*sqrt(w))/10**1.5
       s=ys*sqrt(pi*w*(1-t))/2
       sent=pmax/(b*w*(1-t))
       ax=delk(I)/(1-r)
       else
       if (key.EQ.'MT'. or. key. eq.'mt' )then
       t=2.*ar/w
       sec=1/(cos(pi*t/2))
       ft=sgrt((pi*t*sec)/2)
       snet=pmax/(b*w*(1-t))
       delk(i)=(ft*pp)/(B*sqrt(w))/10**1.5
       ax=delk(I)/(1-r)/1.25
       end if
       end if
       if(ax.ge.s) write(16,98)qq,Ni(gq),Ai(qq),ar,r2,delk(i),dadn(i)
       if(snet.ge.ys) write(16,98)qq,Ni(gq),Ai(qq),ar,r2,delk(i),dadn(i)
       write(16,92)qq,Ni(qq),Ai(qq),ar,r2,dalk(i),dadn(i)
       continue
       write(16,98)qq,Ni(gq),Ai(qq),ar,r2,delk(i),dadn(i)
       continue
       format(i4.lx,i8,2x.f7.3,2x,f7.4,2x,f7.5,2x,f8.4,2x,e10.4)
       format(i4,lx,i8,2x,f7.3,2x,f7.4,2x,f7.5,2x,f7.4,2x,e10.4,2x,'***')
       j=m+4
       k=m+6
       write(16,95)(i,Ni(i),Aai(i),i=j,k)
       format (I4,IX,I8,2X,F7.3)
       write(17,200)(Ni(i),Ai(i),i=l,m+6)
       format(1x,I8,5X,F6.3)
       write (18,350)(delk(i),dadn(i),i=l,qq-3)
       format(1x,F8.5,5X,E12.5)
       close(17)
       close(18)
       close(16)
       close(15)
       end