double precision function proj_test(xi,phi) implicit none integer :: ijk double precision :: xi,phi,xz,yz,rad,cosin,sinus,alpha, & a,b,c,x1,x2,y1,y2,xl,xr,yu,yo, & yl,yr,xo,xu,pi,root_arg,proj,height parameter ( pi = 3.141592653589793d0) ! Test: Kreis mit Mittelpunkt (xz,yz) und Radius rad: xz=0.5 yz=0.5 rad=0.25 if(phi .eq. pi)phi=0.99*pi if(phi .eq. pi/2.d0)phi=0.99*pi/2.d0 cosin=cos(phi) sinus=sin(phi) alpha=xi*sinus + xi*cosin**2/sinus - yz a=1.d0+(cosin/sinus)**2 b=-2.d0*(xz+alpha*cosin/sinus) c=xz**2+alpha**2-rad**2 root_arg=b**2-4.d0*a*c if(root_arg .ge. 0.d0)then x1=(-b-sqrt(root_arg))/(2*a) x2=(-b+sqrt(root_arg))/(2*a) y1=xi*sinus - cosin/sinus*(x1-xi*cosin) y2=xi*sinus - cosin/sinus*(x2-xi*cosin) proj_test=sqrt((x1-x2)**2+(y1-y2)**2) else proj_test=0.d0 endif ! L-Stueck, aus zwei Balken bestehend: do ijk=1,2 if(ijk .eq. 1)then xl=-0.6 xr=-0.4 yu=-0.5 yo=0.7 height=0.75 else xl=-0.4 xr=0.3 yu=-0.5 yo=-0.3 height=0.5 endif cosin=cos(phi) sinus=sin(phi) yl=xi/sinus - cosin/sinus * xl yr=xi/sinus - cosin/sinus * xr ! Es gibt sieben verschiedene Schnittpunkt-Formen: if(yl.gt.yo .and. yr.gt.yo)then proj=0.d0 else if(yl.gt.yo .and. yr.le.yo .and. yr.ge.yu)then xo=sinus/cosin*(xi/sinus-yo) proj=sqrt((xo-xr)**2+(yo-yr)**2) else if(yl.le.yo .and. yl.ge.yu .and. yr.le.yo .and. yr.ge.yu)then proj=sqrt((xr-xl)**2+(yr-yl)**2) else if(yl.gt.yo .and. yr.lt.yu)then xo=sinus/cosin*(xi/sinus-yo) xu=sinus/cosin*(xi/sinus-yu) proj=sqrt((xo-xu)**2+(yo-yu)**2) else if(yl.lt.yu .and. yr.gt.yo)then xo=sinus/cosin*(xi/sinus-yo) xu=sinus/cosin*(xi/sinus-yu) proj=sqrt((xo-xu)**2+(yo-yu)**2) else if(yr.le.yu .and. yl.le.yo .and. yl.ge.yu)then xu=sinus/cosin*(xi/sinus-yu) proj=sqrt((xl-xu)**2+(yl-yu)**2) else if(yl.lt.yu .and. yr.lt.yu)then proj=0.d0 else if(yr.gt.yo .and. yl.le.yo .and. yl.ge.yu)then xo=sinus/cosin*(xi/sinus-yo) proj=sqrt((xl-xo)**2+(yl-yo)**2) else if(yl.lt.yu .and. yr.le.yo .and. yr.ge.yu)then xu=sinus/cosin*(xi/sinus-yu) proj=sqrt((xu-xr)**2+(yu-yr)**2) else print *,'xl xr yo yu yl yr: ',xl,xr,yo,yu,yl,yr stop 'Darf nicht vorkommen!' endif proj_test=proj_test+height*proj enddo end function proj_test