NavList:

   

Reply

The program:

    program cockedhat

c   implicit none

    integer n,i,iseed,itest,in,ratio

    real*4 x(1),sd,h45,h165,h285,pi,z,d2r

    real*4 a,b,c,d,e,g,xx,y,dist,sumsq,sum

    real*4 xsum,xsumsq,ysum,ysumsq,len

    real*4 rratio

    print*,'cockedhat.f'

    pi=4.*atan(1.)

cp  print*,'pi= ',pi

clux

    in=0

    d2r=pi/180.

    ratio=0

    sum=0.

    sumsq=0.

    xsum=0.

    xsumsq=0.

    ysum=0.

    ysumsq=0.

    sd=.5

    print*,'input sd=',sd

    iseed=305

    print*,'iseed= ',iseed

    n=1

clux    type*,'input ii, number of monte-carlo runs'

    print*,'input ii, number of monte-carlo runs'

clux    accept*,ii

    read(5,*)ii

c   ii=10

    print*,'ii',ii

    do 10 j=1,ii

c 45 deg

c y=-tan(45) * x + h45/cos(45)

c y=-1.0 * x + 1.414 * h45

    call NORRAN(N,ISEED,X)

    h45=sd*x(1)

cp  z=cos(pi/180.*45.)

cp  print*,' 45deg,  h45= , y= , x= ',h45,+h45/z,

cp     .  +h45/z/tan(pi/180.*45)



c 165 deg

c y=tan(180-165) * x + -h165/cos(180-165)

c y=.2679 * x - 1.035 * h165

    call norran(n,iseed,x)

    h165=sd*x(1)

cp  z=cos(pi/180.*(180.-165.))

cp  print*,'165deg, h165= , y= , x= ',h165,-h165/z,

cp     .  +h165/z/tan(pi/180.*(180.-165.))



c 285 deg

c y=tan(360-285) * x + h285/cos(360-285)

c y=3.732 * x - 1.035 * h285

    call norran(n,iseed,x)

    h285=sd*x(1)

cp  z=cos(pi/180.*(360.-285.))

cp  print*,'285deg, h285= , y= , x= ',h285,+h285/z,

cp     .  -h285/z/tan(pi/180.*(360-285))



    itest=abs(h45/abs(h45)+h165/abs(h165)+h285/abs(h285))

cc  print*,'itest',itest

    if(itest.EQ.3)in=in+1

cc  print*,'in',in

    a=0.

    b=0.

    c=0.

    d=0.

    e=0.

    g=0.

    a=cos(d2r*45.)**2+cos(d2r*165.)**2+cos(d2r*285.)**2

    b=sin(d2r*45.)*cos(d2r*45.)+sin(d2r*165.)*cos(d2r*165)+

     .      sin(d2r*285.)*cos(d2r*285.)

    c=sin(d2r*45.)**2+sin(d2r*165.)**2+sin(d2r*285.)**2

    d=h45*cos(d2r*45.)+h165*cos(d2r*165.)+h285*cos(d2r*285.)

    e=h45*sin(d2r*45.)+h165*sin(d2r*165.)+h285*sin(d2r*285.)

    g=a*c-b*b

    xx=(a*e-b*d)/g

    y=(c*d-b*e)/g

    dist=sqrt(xx**2+y**2)

cc  print*,j

cc  print*,'xx,y,dist=',xx,y,dist

    sum=sum+dist

    sumsq=sumsq+dist**2

    xsum=xsum+xx

    xsumsq=xsumsq+xx*xx

    ysum=ysum+y

    ysumsq=ysumsq+y*y

    call check(h45,h165,h285,len)

    write(1,900)dist,len,itest

cp  print900,dist,len

900 format(2f12.6,i4)

clux    if(dist .gt. .5887*len)ratio=ratio+1

    rratio=.37135

    if(dist .lt. rratio*len)ratio=ratio+1

10  continue

    print*,' number less than ',rratio,' *len away ratio= ',ratio

    print*,'in, inside cocked hat= ',in

    sumsq=sumsq/ii

    sum=sum/ii

    xsum=xsum/ii

    xsumsq=xsumsq/ii

    ysum=ysum/ii

    ysumsq=ysumsq/ii

    print*,'ii,sumsq,sum',ii,sumsq,sum

    print*,'sd',sqrt(sumsq-sum**2)

    print*,' sum, sumsq,sqrt( sumsq- sum**2)',

     .                 sum,sumsq,sqrt(sumsq-sum**2)

    print*,'xsum,xsumsq,sqrt(xsumsq-xsum**2)',

     .                 xsum,xsumsq,sqrt(xsumsq-xsum**2)

    print*,'ysum,ysumsq,sqrt(ysumsq-ysum**2)',

     .                 ysum,ysumsq,sqrt(ysumsq-ysum**2)

    end



    subroutine check(h45,h165,h285,len)

    real*4 h(3),phi(3),h45,h165,h285,m(3),b(3)

    real*4 x(2),y(2),len

    integer i,ifirst

    phi(1)=45.

    phi(2)=165.

    phi(3)=285.

    h(1)=h45

    h(2)=h165

    h(3)=h285

cp  print*,h

    pi=4*atan(1.)

cp  print*,pi

    do 1 i=1,3

cp  print*,phi(i),-(phi(i)-90.)

    phi(i)=-(phi(i)-90.)*pi/180.

    m(i)=-1./tan(phi(i))

    b(i)=h(i)/sin(phi(i))

cp  print*,i,m(i),b(i)

cp  print'(a3,f6.3,a7,f6.3)','y= ',-1./tan(phi(i)),

cp     .   ' * x + ',h(i)/sin(phi(i))

1   continue

c intersection 1

    x(1)=(b(1)-b(2))/(m(2)-m(1))

    y(1)=x(1)*m(1)+b(1)

cp  print*,x(1),y(1)

    x(2)=(b(1)-b(3))/(m(3)-m(1))

    y(2)=x(2)*m(1)+b(1)

cp  print*,x(2),y(2)

    len=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2)

cp  print*,'len',len

    return

    end



--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to NavList@navlist.net
To , send email to NavList-@navlist.net
-~----------~----~----~----~------~----~------~--~---

   

Reply