174 lines
5.3 KiB
FortranFixed
174 lines
5.3 KiB
FortranFixed
subroutine curvecrossing(ncurves,nsamp,maxnsamp,xdata,ydata,
|
|
¶ms,xcross,ycross)
|
|
implicit none
|
|
!fit ncurves curves from ncurves pairs of datasets and estimate the mean crossings of the curves
|
|
integer ncurves,nsamp(ncurves),maxnsamp,i,j,k,n,m,ndim,
|
|
&iderivative,INFO
|
|
double precision xdata(ncurves,maxnsamp),ydata(ncurves,maxnsamp),
|
|
&xcross,ycross,x(maxnsamp),y(maxnsamp),a(ncurves),b(ncurves),
|
|
&c(ncurves),d(ncurves),beta(10),betamin(10),betamax(10),
|
|
&weitx(maxnsamp),weity(maxnsamp),xmin(maxnsamp),xmax(maxnsamp),
|
|
&y_pred(maxnsamp),x_pred(maxnsamp),sumsquare,
|
|
¶ms(ncurves,4)
|
|
|
|
k=0
|
|
do i=1,ncurves
|
|
do j=1,4
|
|
params(i,j)=-9999.0d0
|
|
enddo
|
|
n=0
|
|
do j=1,nsamp(i)
|
|
if(dabs(xdata(i,j)+9999.0d0).gt.1.0d-8.and.
|
|
&dabs(ydata(i,j)+9999.0d0).gt.1.0d-8)then
|
|
n=n+1
|
|
x(n)=xdata(i,j)
|
|
y(n)=ydata(i,j)
|
|
endif
|
|
enddo
|
|
if(n.gt.1)then
|
|
k=k+1
|
|
beta(1)=30.0d0
|
|
betamin(1)=0.0d0
|
|
betamax(1)=200.0d0
|
|
beta(2)=-4.0d0
|
|
betamin(2)=-100.0d0
|
|
betamax(2)=0.0d0
|
|
beta(3)=8.0d0
|
|
betamin(3)=0.0d0
|
|
betamax(3)=200.0d0
|
|
beta(4)=-0.5d0
|
|
betamin(4)=-10.0d0
|
|
betamax(4)=0.0d0
|
|
do j=1,n
|
|
weitx(j)=1.0d0
|
|
weity(j)=1.0d0
|
|
xmin(j)=0.0d0
|
|
xmax(j)=x(j)+20.0d0
|
|
enddo
|
|
ndim=4
|
|
iderivative=1
|
|
call GenericRegres(n,1,y,1,x,weity,weitx,ndim,beta,betamin,
|
|
&betamax,xmin,xmax,iderivative,0,y_pred,x_pred,sumsquare)
|
|
a(k)=beta(1)
|
|
b(k)=beta(2)
|
|
c(k)=beta(3)
|
|
d(k)=beta(4)
|
|
do j=1,ndim
|
|
params(i,j)=beta(j)
|
|
enddo
|
|
endif
|
|
enddo
|
|
if(k.gt.0)then
|
|
call curmeancrossing(k,a,b,c,d,xcross,ycross)
|
|
else
|
|
xcross=-9999.0d0
|
|
ycross=-9999.0d0
|
|
endif
|
|
return
|
|
end
|
|
|
|
subroutine curmeancrossing(ncurves,a,b,c,d,xcross,ycross)
|
|
implicit none
|
|
!calculate the average crossing point of the ncurves curves.
|
|
integer ncurves,i,j,k
|
|
double precision a(ncurves),b(ncurves),c(ncurves),d(ncurves),
|
|
&xcross,ycross,x(ncurves),y(ncurves),a1,b1,c1,d1,a2,b2,c2,d2,
|
|
&p,q,u,root1,root2,small,term,x1_root1,x2_root1,x1_root2,x2_root2,
|
|
&x_root1,x_root2
|
|
parameter(small=1.0d-7)
|
|
k=0
|
|
do i=1,ncurves-1
|
|
do j=i+1,ncurves
|
|
a1=a(i)
|
|
b1=b(i)
|
|
c1=c(i)
|
|
d1=d(i)
|
|
a2=a(j)
|
|
b2=b(j)
|
|
c2=c(j)
|
|
d2=d(j)
|
|
term=dabs(a1-a2)+dabs(b1-b2)+dabs(c1-c2)
|
|
if(term.lt.small)goto 100
|
|
p=a1-a2+d1-d2
|
|
q=a1*(c2+b1)-a2*(c1+b2)+(d1-d2)*(c2+c1)
|
|
u=a1*b1*c2-a2*b2*c1+(d1-d2)*c1*c2
|
|
call quadraticroots(p,q,u,root1,root2)
|
|
if(root1.gt.0.0d0.or.root2.gt.0.0d0)then
|
|
k=k+1
|
|
if(root1.gt.0.0d0.and.root2.gt.0.0d0)then
|
|
x(k)=root1
|
|
else
|
|
if(root1.gt.0.0d0)then
|
|
x(k)=root1
|
|
else
|
|
x(k)=root2
|
|
endif
|
|
endif
|
|
y(k)=a1*(x(k)+b1)/(x(k)+c1)+d1
|
|
else
|
|
if(dabs(b1-c1).gt.small)then
|
|
term=a2*(b2-c2)/(a1*(b1-c1))
|
|
if(term.ge.0.0d0)then
|
|
term=dsqrt(term)
|
|
root1=((d1+a1)*term+d2+a2)/(term+1.0d0)
|
|
if(dabs(term-1.0d0).gt.small)then
|
|
root2=((d1+a1)*term-d2-a2)/(term-1.0d0)
|
|
else
|
|
root2=-9999.0d0
|
|
endif
|
|
x1_root1=-c1+a1*(b1-c1)/(root1-d1-a1)
|
|
x2_root1=-c2+a2*(b2-c2)/(root1-d2-a2)
|
|
if(dabs(root2+9999.0d0).gt.small)then
|
|
x1_root2=-c1+a1*(b1-c1)/(root2-d1-a1)
|
|
x2_root2=-c2+a2*(b2-c2)/(root2-d2-a2)
|
|
else
|
|
x1_root2=-9999.0d0
|
|
x2_root2=-9999.0d0
|
|
endif
|
|
if(x1_root1.gt.0.0d0.and.x2_root1.gt.0.0d0)then
|
|
x_root1=(x1_root1+x2_root1)/2.0d0
|
|
else
|
|
x_root1=-9999.0d0
|
|
endif
|
|
if(x1_root2.gt.0.0d0.and.x2_root2.gt.0.0d0)then
|
|
x_root2=(x1_root2+x2_root2)/2.0d0
|
|
else
|
|
x_root2=-9999.0d0
|
|
endif
|
|
if(x_root1.gt.0.0d0.or.x_root2.gt.0.0d0)then
|
|
k=k+1
|
|
if(x_root1.gt.0.0d0.and.x_root2.gt.0.0d0)then
|
|
x(k)=x_root2
|
|
y(k)=root2
|
|
else
|
|
if(x_root1.gt.0.0d0)then
|
|
x(k)=x_root1
|
|
y(k)=root1
|
|
else
|
|
x(k)=x_root2
|
|
y(k)=root2
|
|
endif
|
|
endif
|
|
endif
|
|
endif
|
|
endif
|
|
endif
|
|
100 continue
|
|
enddo
|
|
enddo
|
|
if(k.gt.0)then
|
|
xcross=0.0d0
|
|
ycross=0.0d0
|
|
do i=1,k
|
|
xcross=xcross+x(i)
|
|
ycross=ycross+y(i)
|
|
enddo
|
|
xcross=xcross/dble(k)
|
|
ycross=ycross/dble(k)
|
|
else
|
|
xcross=-9999.0d0
|
|
ycross=-9999.0d0
|
|
endif
|
|
return
|
|
end
|