Files
piscal/dataassim/math/othersupmath/linecrossing.f
2016-02-03 18:52:05 +00:00

72 lines
2.0 KiB
FortranFixed

subroutine linecrossing(nlines,nsamp,maxnsamp,xdata,ydata,
&params,xcross,ycross)
implicit none
!fit nlines lines from nlines pairs of datasets and estimate the mean crossings of the lines
integer nlines,nsamp(nlines),maxnsamp,i,j,k,n,m
double precision xdata(nlines,maxnsamp),ydata(nlines,maxnsamp),
&xcross,ycross,x(maxnsamp),y(maxnsamp),intercepts(nlines),
&slopes(nlines),params(nlines,2)
k=0
do i=1,nlines
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
call y_aPLUSbx(n,x,y,intercepts(k),slopes(k))
params(i,1)=slopes(k)
params(i,2)=intercepts(k)
endif
enddo
if(k.gt.0)then
call meancrossing(k,slopes,intercepts,xcross,ycross)
else
xcross=-9999.0d0
ycross=-9999.0d0
endif
return
end
subroutine meancrossing(nlines,slopes,intercepts,xcross,ycross)
implicit none
!calculate the average crossing point of the nlines lines. Parrell lines are excluded
integer nlines,i,j,k
double precision slopes(nlines),intercepts(nlines),xcross,ycross,
&x(nlines),y(nlines),a1,b1,a2,b2
k=0
do i=1,nlines-1
do j=i+1,nlines
a1=slopes(i)
b1=intercepts(i)
a2=slopes(j)
b2=intercepts(j)
if(dabs(a1-a2).gt.1.0d-8)then
k=k+1
x(k)=(b2-b1)/(a1-a2)
y(k)=a1*x(k)+b1
endif
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