C saedigcurve_v2.1.f : reads sae standard fatigue file and creates digital C strain-stress-life curves also in sae std fatigue file format. Jan2000,Mar2005,Mar2009 c Usage: saedigcurve other_filename c linux compile: g77 -g saedigcurve_v2.1.f -o saedigcurve c (in newer Linux compiler comand may be named "gfortran") c Sun Compile: f77 -g -Bstatic saedigcurve_v2.1.f -o saedigcurve C--------------------------------------------------------------------------- C Copyright (C) 1999 SAE Fatigue Design and Evaluation Committee C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of the C license, or (at your option) any later version. C C This program is distributed in the hope tha it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General PUblic License for more details. C C You should have received a copy of the GNU General PUblic License C along with this program; if not, write to the Free Software C Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA C Try also their web site: http://www.gnu.org/copyleft/gpl.html C--------------------------------------------------------------------------- CBugfix Mar4-2009 b In g77 char variables are not pre-set to blanks. Specifically: C the variable used to read #runout Must preset jnp300=" " C Other compilers do pereset to blank. Beware! FAC. Mar4-2009 CBugfix Mar4-2009 Allow for #Fn= tags C Compute Fitted Stress using XKP and XNP rather than Calculated. See CBugFix Mar4-2009 C Remove output of html code (will be done by the calling script "catgnume")facFeb 2009 C Bugfix March 22 2005, MPA version not working correctly. Fix near statm 4007 c Prog to read the SAE standard format, calculate the usual strain-life C parameters and to put out a "fitted" or digital curve file in SAE std format. C Eg. of SAE standard form: C ______ first column in file C | C v C C |# SAE Exchange File Format. C |# data collected from ASTM E606 axial fatigue test data. C | C |#NAME= SAE1045 C |#NAME= SAE350X C |#NAME= SAE050X C |#UNITS= KSI C |#Su= 89. C |#Sy= 50. C |#E= 30000. C |#%RA= 85. C |#BHN= 325 C |#WebPage= http://fde.uwaterloo.ca/Fde/Materials/Steels/ASTM-A588C/g40.21-50A_non_os.html C | C |#note: The data below is not a real file C |# AISI Ph 1, 1998, tested at Waterloo, (iter.1) C |# 0.41%C 1.37Mn .009P .098S .17si .1Cu .05Ni .07Cr .01Mo .002V .001Ca C |# .0029Bo <.011Al .001Ti .0029O .002Cb C |# Rb=91.5 , Upper/Lower Yield points on monotonic C |#LONGITUDINAL to direction of rolling C | C |# Total Strain 2Nf Stress Mean Plastic Strain Initial C |# Amp Amp Stress Amp Elastic Mod. C |0.0125 180 279. .0 0.0030 30100. #specimen comment C |0.0095 490 253. .0 0.0011 29400. C |0.0090 950 229. .0 0.0007 29800. C |0.0075 2260 220. .0 0.0002 30050. C |0.0050 38000. 149. .0 0.0 29900. C |0.0040 770000 119. .0 0 30700. C |# CHARACTER*1 INP1(80) CHARACTER*5 INP5(16),JNP5(16) CHARACTER*10 INP10(8),JNP10(56) CHARACTER*80 JNP80(7),INP80, INP80temp, comment80(50) EQUIVALENCE (INP1(1),INP5(1),INP10(1),INP80), & (JNP5(1),JNP10(1),JNP80(1)) character*300 inp300,jnp300,Cwebpage character*1 inpone(300),jnpone(300) character*10 inpten(30) equivalence (inpone(1), inpten(1), inp300) equivalence (jnpone(1),jnp300) real InMod, longlife integer numcomchars(50) integer*4 iargc, argc real EF, & DE(2500),DS(2500),DEP(2500),DEL(2500),XLIFE(2500), & XMI(2500),XREV(2500),oldDEP(2500),oldDEL(2500) character*10 name1, name2, stressunits, dummy character*30 names30(10), firstfield, ctail character*30 cfiletype, Cdatatype 101 FORMAT(A80) 102 FORMAT(A10) 103 FORMAT(A5) IDIM=2500 XMPAS=6.894759 stressunits="none" C Check to see if a parameter file was specified C Note !!!!!!!!!!!!! for HP fortran the arg numbers must be changed by -1 C and the iragc probably includes the "saehbook" as an arg. c check for correct no of args first Argc = iargc() If (argc .NE. 0) then WRITE(0,*)" saedigcurve: usage ERROR" write(0,*)"Usage e.g.: saedigcurve localFile" STOP ENDIF ninput=0 ndata=0 numnames=0 ifordnumber= 0 numcomment=0 800 continue c Loop back to here for next input line. C Input lines may either be 1. data lines, 2. #Comment lines or 3. blank C It will be hard to distinguish the real comment from the junk comment. C Leave it up to the user to edit the comment section. read(5,"(a300)",end=980)inp300 ninput=ninput+1 Cdebug write(0,*)" read input line ",ninput C Check for blank line if(inp300.eq." ")then C write(6,"(a1)")" " go to 800 endif if(inpone(1).ne."#")then C We may have a data line, or someone screwed up and put the # later in line. C See if 1st char is a # later in line do 805 i=1,300 if(inpone(i).eq." ") go to 805 C No? 1st non-blank found, if # its not a data line loc=i if(inpone(i).eq."#") then C Shift the whole mess over to begining of field inew=0 do 803 j=loc,300 inew=inew+1 inpone(inew)=inpone(j) 803 continue Cdebug write(0,*)"Shifted a comment: ",(inpone(j),j=1,inew) C Ok, now its a nice comment now. Go play with it. go to 880 else C The first non blank is not a # go to 810 endif 805 continue 810 continue C 1st non blank is not a #. Check if its a number. if(inpone(loc).ne."+" .and. & inpone(loc).ne."-" .and. & inpone(loc).ne."." .and. & inpone(loc).ne."0" .and. & inpone(loc).ne."1" .and. & inpone(loc).ne."2" .and. & inpone(loc).ne."3" .and. & inpone(loc).ne."4" .and. & inpone(loc).ne."5" .and. & inpone(loc).ne."6" .and. & inpone(loc).ne."7" .and. & inpone(loc).ne."8" .and. & inpone(loc).ne."9" )then C It must be a letter, not a number. Time to bomb out. write(0,*)" ERROR input line no. ",ninput, & " not a # and not a number. Edit input file" stop endif C Ok, it should be like this one: C 0.0125 180 279. .0 0.0030 30100. #specimen comment C # Total Strain 2Nf Stress Mean Plastic Strain Initial C # Amp Amp Stress Amp Elastic Mod. C Note that the trailing comment field "#specimen comment" may or may not be there. C We should be able to read the first 5 fields as real numbers. Since 2Nf may have C more than 7 signif. digits, it needs to be read in as a double precision. It may C be an integer. The value can be changed later to whatever the local format is. read(inp300,*)StrainAmp, longlife, StressAmp, Smean, & PStrainAmp, InMod Cdebug write(6,*)StrainAmp, longlife, StressAmp, Smean, Cdebug & PStrainAmp, InMod ndata=ndata+1 if(ndata.gt.idim)then write(0,*)" Error too many data points, recompile routine" stop endif C If we have not crashed by here, the data from the line has been read in. C Now figure out if there is comment at the end of line. C Brute force it. Hunt for a # do 820 i=1,300 if(inpone(i).eq."#")then C we found it in col i Cdebug write(0,*)"found # trailer in data line" loc=i go to 822 endif 820 continue C If we got to here, the trailing field is empty go to 855 822 continue C Now count backwards and see where the last char is do 830 i=1,300 j=300-(i-1) if(inpone(j).ne." ")then C found last char lastloc=j go to 832 endif 830 continue 832 continue C It is possible that the #field is an html tag, so it may be real long C Move the stuff to begin of a new field for decoding. Make sure jnp--- is empty CBugfix Mar4-2009 b In g77 char variables are not pre-set to blanks. Make sure here jnp300=" " j=0 do 840 i=loc,lastloc j=j+1 jnpone(j)=inpone(i) 840 continue Cdebug write(0,*)"Trailer # on data: ",(jnpone(i),i=1,j) C Check to see if any special tags are in this field C If the field is blank, then what? read(jnp300,*)ctail Cdebug write(6,*)"#DEBUG: read ",ctail if(ctail.eq."#Runout" .or. & ctail.eq."#RUNOUT" .or. & ctail.eq."#runout" )then C The data point was a runout C In the local companys format, runouts are -ve nos. C Change it to whatever you like !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Cdebug write(6,849)longlife Cdebug 849 format("#FoundRunout=",f20.1) C longlife is a real no. longlife=-longlife endif C #knifeEdge and #outsideGage are also known tags C but we dont do anything with them right now. C So just write out, and say which point. Cdebug write(6,850)ndata,(jnpone(i),i=1,j) 850 format(" Pt No.", I5, " : ", 300a1) C Now save the data in the local program fields 855 continue DE(ndata)=StrainAmp DS(ndata)=StressAmp DEP(ndata)=PStrainAmp DEL(ndata)=StrainAmp-PStrainAmp XLIFE(ndata)=longlife XMI(ndata)=StressAmp/DEL(ndata) XREV(ndata)=0. C oldDEP and oldDEL are not used right now. C End of data line processing, write it out later, when all are in. go to 800 endif C Input line has a # in first col, check it just in case 880 continue if(inpone(1).ne."#")then C something is bad in program write(0,*)" ERROR 880, sorry prog. messed up call ? admin" stop endif C Ok, its a nice comment. Figure out if its a special tag. read(inp300,*)firstfield if(firstfield .eq."#FileType=" .or. & firstfield .eq."#filetype=" .or. & firstfield .eq."#FILETYPE=" )then read(inp300,*) firstfield, cfiletype if(cfiletype.ne."strain_life")then C We have the wrong kind of sae file here folks write(0,*)" ERROR, wrong type of SAE std file. ", & " Not #FileType= strain_life" stop endif C filetype is ok. Write it out as a comment go to 950 endif if(firstfield .eq."#DataType=" .or. & firstfield .eq."#DATATYPE=" .or. & firstfield .eq."#datatype=" )then read(inp300,*) firstfield, Cdatatype if(Cdatatype.ne."raw")then C We have the wrong kind of sae file here folks write(0,*)" ERROR, wrong type of SAE std file. ", & " Not #DataType= raw" stop endif C DataType is ok. Write it out as a comment C Change it to "fitted", since that is what we are doing in C this progr. write(inp300,888) 888 format("#DataType= fitted") go to 950 endif if(firstfield .eq."#NAME=" .or. & firstfield .eq."#Name=" .or. & firstfield .eq."#name=" )then numnames=numnames+1 read(inp300,*) firstfield, names30(numnames) go to 950 endif if(firstfield .eq."#UNITS=" .or. & firstfield .eq."#units=" .or. & firstfield .eq."#Units=" .or. & firstfield .eq."#Stress_units=" .or. & firstfield .eq."#STRESS_UNITS=" .or. & firstfield .eq."#Stress_Units=" .or. & firstfield .eq."#stress_units=" )then read(inp300,*) firstfield, stressunits go to 950 endif if(firstfield .eq."#Su=" .or. & firstfield .eq."#SU=" )then read(inp300,*) firstfield, Sult go to 950 endif if(firstfield .eq."#Sy=" .or. & firstfield .eq."#SY=" )then read(inp300,*) firstfield, Syield go to 950 endif if(firstfield .eq."#E=" .or. firstfield .eq."#e=" .or. & firstfield.eq."#EMOD=" .or. firstfield .eq."#Emod=".or. & firstfield .eq."#emod=" .or. & firstfield .eq."#MODULUS=" .or. & firstfield .eq."#Modulus=" .or. & firstfield .eq."#modulus=" )then read(inp300,*) firstfield, EMOD go to 950 endif if(firstfield .eq."#%RA=" .or. & firstfield .eq."#%Ra=" )then read(inp300,*) firstfield, percentRA go to 950 endif if(firstfield .eq."#WebPage=" .or. & firstfield .eq."#Webpage=" .or. & firstfield .eq."#WEBPAGE=" .or. & firstfield .eq."#webpage=" )then read(inp300,*) firstfield, Cwebpage go to 950 endif if(firstfield .eq."#FractureStress=" .or. & firstfield .eq."#fracturestress=" .or. & firstfield .eq."#FRACTURESTRESS=" .or. & firstfield .eq."#Fracturestress=" )then read(inp300,*) firstfield, FractureStress go to 950 endif if(firstfield .eq."#FractureStrain=" .or. & firstfield .eq."#Fracturestrain=" .or. & firstfield .eq."#fracturestrain=" .or. & firstfield .eq."#FRACTURESTRAIN=" )then read(inp300,*) firstfield, FractureStrain go to 950 endif if(firstfield .eq."#Ford=" .or. & firstfield .eq."#FORD=" .or. & firstfield .eq."#ford=" .or. CBugfix Mar4-2009 Allow for #Fn= & firstfield .eq."#Fn=" .or. & firstfield .eq."#fn=" )then read(inp300,*) firstfield, ifordnumber go to 950 endif if(firstfield .eq."#BHN=" .or. & firstfield .eq."#bhn=" .or. & firstfield .eq."#Bhn=" .or. & firstfield .eq."#HBN=" .or. & firstfield .eq."#HB=" )then read(inp300,*) firstfield, fbrinell nbrinell=ifix(fbrinell) goto 950 endif C Look for things we should skip over in the output write(INP80temp,939) 939 format("#Here is the bottom part of the ", & "html graph/calc wrapper-----------") if(INP300.eq. INP80temp)go to 800 if(INP300.eq. &"#" & )go to 800 if(INP300.eq. &"#" & )go to 800 CC It was NOT a special comment. Save it. If it was a # blank line, discard it. CC Count backwards to find last non-blank char of string C do 940 i=1,300 C j=300-(i-1) C if(inpone(j).ne." ") then CC found last char C lastloc=j C go to 943 C endif C 940 continue C CC We also need to save the no of chars for each comm line CC If it was a #blank forget it C 943 continue C if(lastloc.eq.1)go to 800 CC Look for some other dumb lines and discard them C read(inp300,*) dummy C if(dummy.eq."#saeinput")go to 800 C if(dummy.eq."#---------")go to 800 C C numcomment=numcomment+1 C numcomchars(numcomment)=lastloc CC save the first 80 chars C read(inp300,"(a80)")comment80(numcomment) C go to 800 950 continue C Write out the comment line C Count backwards and see where the last char is do 960 i=1,300 j=300-(i-1) if(inpone(j).ne." ")then C found last char lastloc=j go to 962 endif 960 continue 962 continue write(6,"(300a1)")(inpone(i),i=1,lastloc) C Go read another line go to 800 980 continue C All input lines have been read. The comments were C put out along the way. Its time to dump out the data C in whatever format the local machine wants it. This bit C is probably site specific. if(stressunits .eq."mpa" .or. & stressunits .eq."MPA" .or. & stressunits .eq."MPa" .or. & stressunits .eq."Mpa" )then stressunits="MPa" C Save this for later reconversion to MPa, if required. endif if(stressunits .eq."ksi" .or. & stressunits .eq."KSI" .or. & stressunits .eq."Ksi" )then stressunits="KSI" endif EMODmpa=EMOD C This local machine s data runs in ksi, change to ksi if(stressunits .eq. "MPa")then EMOD=EMOD/XMPAS Sult=Sult/XMPAS Syield=Syield/XMPAS FractureStress=FractureStress/XMPAS do 983 i=1,ndata DS(i)=DS(i)/XMPAS XMI(i)=XMI(i)/XMPAS 983 continue endif C Output for a datafile (not for saehandbook (old FD+E work) C WRITE(6,*)BHNnumber,Emod,0.,0.,0.,Syield,0. C WRITE(6,*)0.,Sult,0.,0., 0.,0.,0. C WRITE(6,*)0.,0.,0.,percentRA,0.,FractureStress,0. C WRITE(6,*)FractureStrain C WRITE(6,*)ndata C C DO 988 I=1,ndata C WRITE(6, 978)DE(I),DS(I),DEP(I),DEL(I),XLIFE(I), C & XMI(I),XREV(I) C 978 FORMAT("#"1X,F9.6,1X,F8.2,2(1X,F9.6),1X,F11.0,1X, C & F8.0,1X,F7.2) C 988 continue C analysis & plotting------------------------------------------------ 3770 CONTINUE XMODM=EMOD*XMPAS/1000. YLDM=Syield*XMPAS SUM=Sult*XMPAS XKM=XK*XMPAS STFM=FractureStress*XMPAS NUMF=ndata XMOD=EMOD FLAT=0. C =0 means no "Flat Top Yield", =Syield then assume flat top yield IYNP=+1 C = Plot pts on plastic & elastic vs life C =-1 no. then do not plot points IDIST=0 C IDIST=1 FOR NORMAL DISTRIBUTION OF TOLERANCE BOUNDS C =0 FOR STUDS T DISTB. OF BANDS C Run the data fitting so that we can print the parameters: CALL LEAST(IDIST,DEP,DS,IDIM,NUMF,XKP,XNP,STD, & C1,XKU,XKL,R,T,IERR) CALL LEAST(IDIST,XLIFE,DEP,IDIM,NUMF,EFP,C,STD,C1, & BU,BL,R,T,IERR) CALL LEAST(IDIST,XLIFE,DEL,IDIM,NUMF,SFPE,BED,STD, & C1,BUE,BLE,R,T,IERR) SFP=SFPE*XMOD SFPM=SFP*XMPAS XKPM=XKP*XMPAS Compute the "Corrected values" of K' and n' CORXN=0. CORXK=0. IF(BED.NE.0. .AND. C.NE.0.)CORXN=BED/C IF(CORXN.NE.0. .AND.EFP.NE.0.)CORXK=SFP/(EFP**CORXN) CORXKM=CORXK*XMPAS C Calc cyclic yield strength 0.2% offset, both the experimental value, C and the fitted curve value. CYCYIELD=CORXK*(0.002)**CORXN CYCYIELDM=CYCYIELD*XMPAS CYCYLDX=XKP*(0.002)**XNP CYCYLDXM=CYCYLDX*XMPAS C Write out the title line C Get the first 10 chars from each of first two names name1=" " name2=" " if(numnames.ge.2)then read(names30(1),*)name1 read(names30(2),*)name2 endif if(numnames.eq.1)then read(names30(1),*)name1 endif write(6,1490),name1,name2, nbrinell,ifordnumber 1490 format("# ",A10," ",A10," BHN= ",i4," Fn= ",I4) CC Now write out a bunch Comment lines below the title C do 1500 i=1,numcomment C INP80=comment80(i) C CC Look for things we should skip over in the output C write(INP80temp,1492) C 1492 format("#Here is the bottom part of the ", C & "html graph/calc wrapper-----------") C if(INP80.eq. INP80temp)go to 1500 C C if(INP80.eq. C &"#" C & )go to 1500 C C if(INP80.eq. C &"#" C & )go to 1500 C C write(6,"(80a1)" )(inp1(j),j=1,numcomchars(i)) C 1500 continue C C C 123456789_123456789_123456789_123456789_123456789_123456789_ write(6,1501) 1501 format("# Monotonic Props. Cyclic Props.") write(6,1510)EMOD,XMODM,CORXK,CORXKM 1510 format("#ELAS. MOD.= ",F6.0," KSI, ",F5.0," GPA", &T40,"K'",T50,"=",F6.1," KSI,",F7.0,"MPA") CalKprime=CORXK CalNprime=CORXN write(6,1520)Syield, YLDM, CORXN 1520 format("#YIELD,0.2%= ",F6.0," KSI, ",F5.0," MPA", &T40,"N'",T50,"=",F7.4) write(6,1530)Sult,SUM,SFP,SFPM 1530 format("#ULT. STRG.= ",F6.0," KSI, ",F5.0," MPA", &T40,"F. STRG COEF=",F7.1," KSI,",F6.0,"MPA") write(6,1540)XK,XKM,BED 1540 format("#K = ",F8.1," KSI, ",F7.0," MPA", &T40,"F.STRG EXP, b=",F7.4) write(6,1550)XN,EFP 1550 format("#N",12X,"=",F7.4, T40,"FAT DUCT COEF=",F9.4) write(6,1560)percentRA,C 1560 format("#RED. IN AREA =",F5.1,T40,"F.DUCT EXP, c=",F7.4) SFracTrue=0. write(6,1570)SFracTrue,STFM,CYCYIELD,CYCYIELDM 1570 format("#T. FRAC. STG.= ",F5.1," KSI, ",F5.0," MPA", &T40,"Exp Cyc Yld =",F6.0," Ksi,",F5.0,"MPA") EF=0. write(6,1580)EF,CYCYLDX,CYCYLDXM 1580 format("#T. FRAC. STR.= ",F5.3, &T40,"Fit Cyc Yld =",F6.0," Ksi,",F5.0,"MPA") write(6,1590)ndata 1590 format("#No. fatigue data points=", I5) C ------------------------------------ C PRINT out the "fitted" strain stress life curve data C C eg.format: C |# Total Strain 2Nf Stress Mean Plastic Strain Initial C |# Amp Amp Stress Amp Elastic Mod. C |0.0125 180 279. .0 0.0030 30100. #fitted_point C write(6,4004) 4004 format("#") write(6,4005) 4005 format("# NOTE!! The Following Points are ", & "FITTED DATA:"/ & "#NOTE!! Fitted Stress computed using Experm. K' and n'" &) if(stressunits .eq. "KSI") write(6,4006) 4006 format("#Stress_Units= KSI") if(stressunits .eq. "MPa") write(6,4008) 4008 format("#Stress_Units= MPa") if(stressunits .ne. "MPa" .and. & stressunits .ne. "KSI")then write(6,4009) 4009 format("#NOTE!! No stress units specified! Assuming:") write(6,4006) stressunits = "KSI" endif write(6,4003) 4003 format("# Total Strain 2Nf Stress Mean ", & "Plastic Strain Initial"/ & "# Amp Amp Stress ", & "Amp Elastic Mod." ) SF=SFP/EMOD YTOP=SF+ EFP XTOP=1.0 Smean=0. CBugFix Mar4-2009 Use the "experimental" K' and n' to compute C the Fitted Stress amplitude C CALL getstress(EMOD,CalKprime,CalNprime,YTOP, Stress) CALL getstress(EMOD,XKP,XNP,YTOP, Stress) Eplastic= YTOP-(Stress/EMOD) IXTOP=ifix(XTOP) Eout=EMOD if(stressunits.eq."MPa")then C The original file stress was MPa, change them back here Eout=EMODmpa Stress=Stress*XMPAS endif write(6,4007)YTOP,IXTOP, Stress,Smean, Eplastic,Eout 4007 format(1x, f9.5,1x,i10, 1x,f7.1,1x,f7.0, 1x,f9.5,4x,f7.0, & " #Fitted_point") C Loop from decade to decade and print out some of the curve s values C The procedure here is to get a life -> strain -> stress 4015 continue C DRAW TO 1/3 DECADE XN2=XTOP*2.0 C Enough? (10**8) IF(XN2 .GT. 1.0E9)go to 4020 C Nope, more yet YN2=SF*(XN2)**BED + EFP*(XN2)**C CBugFix Mar4-2009 Use the "experimental" K' and n' to compute C the Fitted Stress amplitude C CALL getstress(EMOD,CalKprime,CalNprime,YN2, Stress) CALL getstress(EMOD,XKP,XNP,YN2, Stress) Eplastic= YN2-(Stress/EMOD) IXTOP=ifix(XN2) if(stressunits.eq."MPa")Stress=Stress*XMPAS write(6,4007)YN2,IXTOP, Stress,Smean, Eplastic,Eout C Now go to 2/3 of Decade XN3=XTOP*5.0 YN3=SF*(XN3)**BED + EFP*(XN3)**C CBugFix Mar4-2009 Use the "experimental" K' and n' to compute C the Fitted Stress amplitude C CALL getstress(EMOD,CalKprime,CalNprime,YN3, Stress) CALL getstress(EMOD,XKP,XNP,YN3, Stress) Eplastic= YN3-(Stress/EMOD) IXTOP=ifix(XN3) if(stressunits.eq."MPa")Stress=Stress*XMPAS write(6,4007)YN3,IXTOP, Stress,Smean, Eplastic,Eout XNEXT=XTOP*10.0 C See if we have enough IF (XNEXT .GT. 1.0E9) go to 4020 C not yet YNEXT= SF*(XNEXT)**BED +EFP*(XNEXT)**C CBugFix Mar4-2009 Use the "experimental" K' and n' to compute C the Fitted Stress amplitude C CALL getstress(EMOD,CalKprime,CalNprime,YNEXT, Stress) CALL getstress(EMOD,XKP,XNP,YNEXT, Stress) Eplastic= YNEXT-(Stress/EMOD) IXTOP=ifix(XNEXT) if(stressunits.eq."MPa")Stress=Stress*XMPAS write(6,4007)YNEXT,IXTOP, Stress,Smean, Eplastic,Eout XTOP=XNEXT YTOP=YNEXT GO TO 4015 4020 continue C 9000 CONTINUE STOP END C--------------------------------------------------------- SUBROUTINE getstress(CEMOD,xkprime,xnprime,Strain, Stress) C S/r computes Stress, given Strain, etc. C CEMOD=Elast. Mod, Stress & Strain are amplitudes XINNP=1.0/xnprime XINVE=1.0/CEMOD EXPO=XINNP-1.0 CKP=xkprime E=Strain C INTERPOLATION REQUIRED SOLD=CKP Coct97 A=(1.0/CKP)**XINNP C With some data sets, we get problems here when the Fat. Duct Coeff is very high C On CCChus suggestions we changed the code here, and in the iteration below C Oct.17, 1996 fac DO 6520 I2=1,1000 Coct97 RISE=SOLD*XINVE +A*(SOLD)**(XINNP)-E Coct97 DERIV=XINVE+A*XINNP*SOLD**(EXPO) RISE=SOLD*XINVE + (SOLD/CKP)**(XINNP)-E DERIV=XINVE + XINNP/CKP*(SOLD/CKP)**(EXPO) SNEW=SOLD-(RISE/DERIV) IF (ABS(SOLD-SNEW)/SNEW .LT. 0.00001) GO TO 6525 SOLD=SNEW 6520 CONTINUE 6525 Stress=SNEW return end C========================================== SUBROUTINE LEAST(IDIST,X,Y,IM,N,BE,M,STD,C1,C2E,C2E2, + R,T,IERR) C* C* S/R FOR PERFORMING LEAST SQUARES ON LOG-LOG BASIS. C* ADOPTED FROM LANDGRAF'S BASIC ROUTINE. C* SEPT.7/79 A.C. C* C*FORMULA ASSUMED: Y=MX + B C* OR LOGY = M* LOGX + B C*X,Y = LIFE AND STRAIN VECTORS C*N = NO OF DATA PTS IN X & Y C*IM = DIMENSION OF X&Y VECTORS C*FN = FLOAT PT. OF N LESS RUNOUTS (RUNOUTS -VE X'S) C*B = LOG-LOG INTERCEPT C*BE = EXP(B) , COEFF. OF POWER LAW C*C2 = ERROR ON B C*C2E,C2E2 = UPPER AND LOWER ERROR ON "BE" C*M = SLOPE OF LINE, EXP. OF POWER LAW C*C1 = ERROR ON SLOPE, M C*R = CORR. COEFF. C*T = STUDENT'S T FOR FN SPECS. 95% C*IDIST=1 FOR NORMAL DISTRIBUTION OF TOLERANCE BOUNDS C* =0 FOR STUDS T DISTB. OF BANDS C*IERR = ERROR INDICATOR C* REAL X(IM), Y(IM), XL(500),YL(500),M,XT(10) DATA XT/12.706, 4.3303, 3.182, 2.776, 2.571, 2.447, + 2.365, 2.306, 2.262, 2.228/ C* IERR=0 IF (N .LT. 2) GO TO 20 IF (N .GT. 500) GO TO 30 C* COMPUTE LEAST SQS. C* CHANGE TO LN(X), AND LN(Y) FIRST. SUMX=0.0 SUMY=0.0 NUMB=0 DO 60 I=1,N C* CHECK FOR RUNOUTS, IF -VE , ELIMINATE C* ALSO IF PLASTIC STRAIN (VS. STRESS) IS LESS THAN 0.0001 IF (X(I) .LE. 0.00010) GO TO 60 IF (Y(I) .LE. 0.0001) GO TO 60 NUMB=NUMB+1 XLT=ALOG(X(I)) XL(NUMB)=XLT SUMX=SUMX+XLT YLT=ALOG(Y(I)) YL(NUMB)=YLT SUMY=SUMY+YLT 60 CONTINUE C* Changed Oct.96, to allow 2 pts to be fit Change IF(NUMB.LE.2) GO TO 20 IF(NUMB.LT.2) GO TO 20 FN=FLOAT(NUMB) XBAR=SUMX/FN YBAR=SUMY/FN CC=0.0 DD=0.0 EE=0.0 DO 70 I=1,NUMB XLT=XL(I)-XBAR YLT=YL(I)-YBAR CC=CC+XLT*YLT DD=DD+ XLT**2 EE=EE + YLT ** 2 70 CONTINUE C* C* COMPUTE FACTORS M=CC/DD B= YBAR- M * XBAR BE= EXP(B) C* COMPUTE STANDARD DEVIATION SUM=0.0 DO 80 I=1,NUMB EXPY=B+M*XL(I) 80 SUM=SUM+(EXPY-YL(I))**2 STD=SQRT(SUM/(FN-2)) C* COMPUTE STUDENT'S T FOR 95% CONFID. INTERV. C* NM2= NO OF DEG. OF FREEDOM NM2=NUMB-2 IF (NM2 .GT. 1) GO TO 6 C* NOT ENOUGH POINTS FOR STD'S T IERR=2 T=0.0 GO TO 50 6 CONTINUE C* CHECK FOR NORMAL OR STUD'S T DISTRIBUTION REQ. T=1.96 IF (IDIST .EQ. 1 ) GO TO 50 IF (NM2 .GT. 300) GO TO 5 IF (NM2 .LE. 10) GO TO 1 IF (NM2.LE. 20) GO TO 2 IF (NM2 .LE. 40) GO TO 3 IF (NM2 .LE.100) GO TO 4 C* NM2 .GT. 100 & LT. 300, INERPOLATE T=1.984-(1.984-1.9679)*FLOAT(NM2-100)/200. GO TO 50 4 T=2.021-(2.021-1.984)* FLOAT( NM2-40 )/ 60.0 GO TO 50 3 T= 2.086- (2.086-2.021) *FLOAT(NM2-20)/ 20.0 GO TO 50 2 T= 2.228- (2.228-2.086) *FLOAT(NM2-10)/ 10.0 GO TO 50 1 T=XT(NM2) GO TO 50 5 T=1.9679 GO TO 50 C* C* NOT ENOUGH PTS. 20 CONTINUE C NOT ENOUGH PTS FOR L SQ. B=0. BE=0. M=0. IERR=1 RETURN 30 WRITE (6,31) N 31 FORMAT(" TOO MANY PTS, CHANGE DIMENSION OF VECTORS"/ + " XL & YL IN S/R LEAST. N= ",I6/) IERR=2 RETURN C* 50 CONTINUE C1= T*STD C2=C1 C2E= EXP(B+C2) C2E2= EXP(B-C2) R= CC/(SQRT(DD*EE)) IERR=0 RETURN END C============================================