Hjælp til Fortran
HejEr 100% ny inden for feltet, men skal bruge program som kan køre nedenstående kode. Jeg har forsøgt med Force 2.0.8 men den kender ikke alle routinerne. der er 2 tilhørende filer med. det skal helst være et program til Win.
parameter(nmax=2000)
real grams,sec,time,force,delt,test
dimension grams(nmax),sec(nmax)
character*256 dum1,dum2,dum3,sample,solvent,comment
real break(nmax), cscoef(4,nmax), fdata(nmax), xdata(nmax),x
real smpar,sdev,weight(nmax)
real y(4,nmax)
integer nintv
external csder,cssmh, csval, sset
integer jtest,ntot,jj
* parameters for the plotting routine
character*256 xlabel,ylabel,plotlabel
integer numpts
real plxmin,plxmax,plymin,plymax
real plotx(nmax),ploty(nmax)
open(unit=7,file='raw.dat',status='old')
open(unit=8,file='rectfit.dat',status='unknown')
open(unit=10,file='prop.dat',status='old')
* discard header info
do j=1,6
read(7,*)dum1
end do
* read sample name, solvent, width, depth, deflection and span of plate (mm)
* here, a and b are full thickness and width, rather than half-thicknesses
read(7,*)dum1,dum2,dum3,sample
read(7,*)dum1,dum2,solvent
read(7,*)dum1,dum2,dum3,b
read(7,*)dum1,dum2,dum3,a
read(7,*)dum1,dum2,dum3,delta
read(7,*)dum1,dum2,dum3,span
* write sample name and solvent into output file
write(8,*)sample
write(8,*)solvent
* convert to cm
span = span/10.
delta = delta/10.
a = a/10.
b = b/10.
* discard column titles
read(7,*)dum1
* read elapsed time, load, displacement, jtest, and temperature
* (jtest=5 if load is applied, =0 otherwise)
read(7,*)time,force,delt,test,temp
jtest=int(test)
sec(1)=time
grams(1)=force
* sum is baseline load
sum=grams(1)
* dsum0 is baseline value of displacement (mm)
dsum0=delt
* choose initial values for time and force (before bending bar)
n=2
5 read(7,*)time,force,delt,test,temp
jtest=int(test)
if(jtest.eq.0)then
dsum0 = dsum0+delt
sec(n)=time
grams(n)=force
sum=sum+grams(n)
n=n+1
if(n.ge.nmax)goto 95
go to 5
end if
* gram0 = baseline load
gram0 = sum/(n-1.)
* dsum0 = baseline displacement
dsum0=dsum0/(n-1.)
* n0 = number of points before load is applied
n0 = n
write(6,*)'Number of points in baseline =',n0
* discard line in which jtest becomes > 0
* discard data for time < discard
read(10,*)comment
read(10,*)discard
10 read(7,*,end=95)time,force,delt,test,temp
if(time.lt.discard)go to 10
* dsum is displacement (mm) during loading
dsum=0.
jj=0
joff = 0
20 read(7,*,end=30)time,force,delt,test,temp
jtest=int(test)
n = n+1
if(n.ge.nmax)goto 95
* jtest = 5 during the measurement, then becomes 0 when the load is lifted
if(jtest.ne.0) then
* adjust sec to be nonzero immediately upon application of load
* subtract baseline value from load so grams = applied load in grams
sec(n)=time+0.1
grams(n)=gram0-force
dsum=dsum+delt
jj=jj+1
* write values into arrays to be plotted
plotx(jj)=real(log10(sec(n)))
ploty(jj) = grams(n)
go to 20
else
* if the jtest = 0, the experiment is over.
* there are ntot data sets
ntot = jj
end if
write(6,*)'ntot =',ntot
* discard the first line after test changes
read(7,*,end=30)time,force,delt,test,temp
* read data for time after load applied
sum=0.
joff = 0
* average the final baseline to check for drift of the load cell
30 read(7,*,end=35)time,force,delt,test,temp
sum=sum+force
joff=joff+1
go to 30
35 drift = sum/(joff)-gram0
* subtract baseline from displacement and convert to cm
dsum=(dsum0-dsum/ntot)/10.
if(abs(1.-delta/dsum).gt.0.001) then
write(6,*)' actual delta =',dsum,', nominal delta =',
&delta
delta=dsum
else
continue
end if
* write values into rectfit.dat
write(8,100)span,a,b,delta
write(6,120)drift
nstart = n0+2
nstop = nstart+ntot-2
nintv = 0
do j=nstart,nstop
nintv = nintv + 1
xdata(nintv)=alog10(sec(j))
fdata(nintv)=grams(j)
end do
* smooth the raw data
sdev = 1.e-4*abs(fdata(1))
call sset (nintv, sdev, weight, 1)
smpar = nintv
call cssmh (nintv, xdata, fdata, weight, smpar, break,cscoef)
* find the first derivative of the smoothed data
do j=1,nintv
x = xdata(j)
* y(1,j) contains the smoothed values of y = w(t)
y(1,j)=csval(x,nintv-1,break,cscoef)
* fdata now contains the first derivative, y'
fdata(j) = csder(1,x,nintv,break,cscoef)
end do
* smooth y'
sdev = 1.e-4*abs(fdata(1))
call sset (nintv, sdev, weight, 1)
smpar = nintv
call cssmh (nintv, xdata, fdata, weight, smpar, break,cscoef)
* find the second derivative, y''
do j=1,nintv
x = xdata(j)
* y(2,j) contains the smoothed values of y'
y(2,j)=csval(x,nintv-1,break,cscoef)
* fdata now contains the second derivative, y''
fdata(j) = csder(1,x,nintv,break,cscoef)
end do
* smooth y''
sdev = 1.e-4*abs(fdata(1))
call sset (nintv, sdev, weight, 1)
smpar = nintv
call cssmh (nintv, xdata, fdata, weight, smpar, break,cscoef)
* find the maximum in y'', which corresponds to the inflection point of interest
ymax=0.
do j=1,nintv
x = xdata(j)
* y(3,j) contains the smoothed values of y''
y(3,j)=csval(x,nintv-1,break,cscoef)
if((ymax.gt.0.).and.(y(3,j).lt.0.))go to 40
if(y(3,j).gt.ymax) then
ymax=y(3,j)
tau=x
end if
end do
40 tau=5.*10**tau
write(6,*)' tau =',tau
* delta is approximate loading period
delta = -0.5*(sec(nstart+1)-sec(nstart))
write(6,*)' volume fraction solids, bulk modulus of fluid ',
& '(gpa),'
write(6,*)' and estimates of exponent m and poisson ratio'
read(10,*)comment
read(10,*)rho,bmf,ex,pr
write(6,*)rho,bmf,ex,pr
write(6,*)' Hit carriage return to see plot'
* estimate initial load from first point after application (to allow for bounce)
theta = sec(nstart)/tau
* ratio = w(theta)/w(0), assuming a = 0.05
ratio = 0.9+0.05*exp(-4.51*(theta**0.5-theta**2.5)/
& (1.-theta**0.575))
w0 = grams(nstart)/ratio
write(8,100)delta,w0,tau
write(8,100)rho,bmf,ex,pr
do j=nstart,nstop
write(8,100)sec(j),grams(j)
end do
pause
********************************************************************
* set up plot
call plsdev("Mac1")
call plscolbg(255,255,255)
call plinit()
* specify number of points
numpts = ntot
* specify limits of axes
plxmin = plotx(1)
plxmax = plotx(numpts)
plymin = ploty(numpts)
plymax = ploty(1)
call plcol(9)
call plenv(plxmin,plxmax,plymin,plymax,0,0)
* provide labels for axes and for plot
xlabel = 'log(t) (seconds)'
ylabel = 'W(t) (grams)'
plotlabel = 'Relaxation data (Hit any key to continue)'
call plcol(8)
call pllab(xlabel,ylabel,plotlabel)
* plot points
* plot individual points with symbol corresponding to 6
* set color of points to black
call plcol(8)
call plpoin(numpts,plotx,ploty,1)
* close the plotting routine
call plend()
stop
95 write(6,*)' too many data'
stop
100 format(1x,5(e12.5,5x))
110 format(1x,2(i6,3x))
120 format(5x,'baseline drift (g) =',f9.4)
end
på forhånd tak
