Welcome to the NavList Message Boards.


A Community Devoted to the Preservation and Practice of Celestial Navigation and Other Methods of Traditional Wayfinding

Compose Your Message

Add Images & Files
    Cocked Hat Monte Carlo Program
    From: Dave Walden
    Date: 2007 Mar 17, 08:33 -0700
    Repeating, so it's all together.  I've successfully compiled with g77 under linux and f77 (an old DEC/HP/MS compiler) under NT.
    Below is a short FORTRAN program to do Monte Carlo simulations of a three point LOP.  It is rough and uncommented.  I provide it as is since there seems to be quite some interest at the moment.  I have found it, and other versions and modifications, very helpful as others have said recently such a program might be.  It's late, and I apolgize in advance if I mis-speak below.
    It assumes, without loss of generality, the true location is at x=0, y=0.  sd is the standard deviation of altitude uncertainty.  For each Monte Carlo run, for each LOP it finds the altitude error based on a random selection from a distribution of errors. h45 for example. (NORRAN is the subroutine.  The choice of random number generators can require some attention, as many may realize.)  It checks if the true location is inside the cocked hat, itest=3.  It finds the center of the cocked hat, the a,b,c,d,e,g calculation from the Nautical Almanac.  It calculates the distance from the true location to the center of the cocked hat, dist.  It accumulates some variables to provide statistics on this distance as well as x and y distance from the true location. It calculates the length of the side of this cocked hat, len, via CALL CHECK.  It calculates the number of times the point is with a cirle of radius of .5887*len of the true location.  It writes out dist, len, and itest.  By analyzing this file, one can calculate probabilites of such things as small cocked hats that do not contain the true location, how often large cocked hats do contain cocked hats, how often cocked hats of various sizes occur, etc.
    These quantities can be calculated analytically and I have found the Monte Carlo results to tend towards these results as the number of runs increases.  (A million only takes a minute on an old, slow PC.)
    In the sample output, next, one sees that of the 100,000 runs for 25,230 cases, the true location was inside the cocked hat. Approaching 25%.   The average x and y errors are -.0012 and .0008 respectively.  Approaching zero as expected.
     ii 100000   number less than   0.37134999 *len away ratio=  26119  in, inside cocked hat=  25230  ii,sumsq,sum 100000  0.331236064  0.510192156  sd  0.26634568   sum, sumsq,sqrt( sumsq- sum**2)  0.510192156  0.331236064 0.26634568  xsum,xsumsq,sqrt(xsumsq-xsum**2) -0.00127442868  0.16588138 0.407283396  ysum,ysumsq,sqrt(ysumsq-ysum**2)  0.0008370672  0.165351331 0.406633288 sh-3.00#
    The first few lines from the output file,
        0.057159    1.233867   3
         0.724052    0.217705   1
         0.387911    0.322774   1
         0.608533    0.352897   1
         0.635789    0.097189   1
         1.136409    1.546760   1 
        0.175892    0.142276   1
    show for example, a case where the distance from the center of the cocked hat to the true location is .057 min (or naut miles) and length of the sides of the cocked hat is 1.233 min and the true location is inside the cocked hat, itest=3.

    Don't be flakey. Get Yahoo! Mail for Mobile and
    always stay connected to friends
    To post to this group, send email to NavList@fer3.com
    To unsubscribe, send email to NavList-unsubscribe@fer3.com

    Browse Files

    Drop Files


    What is NavList?

    Join NavList

    (please, no nicknames or handles)
    Do you want to receive all group messages by email?
    Yes No

    You can also join by posting. Your first on-topic post automatically makes you a member.

    Posting Code

    Enter the email address associated with your NavList messages. Your posting code will be emailed to you immediately.

    Email Settings

    Posting Code:

    Custom Index

    Start date: (yyyymm dd)
    End date: (yyyymm dd)

    Visit this site
    Visit this site
    Visit this site
    Visit this site
    Visit this site
    Visit this site