External Function and Program Repository (K-O)

 

| HOME |   BACK   |

LEASTLEV.prg

 

/ ************************************************************************

/                Least Square Leveling Network Adjustment

/ written by : Chainman
/ date       : 12/05/2001

/ ************************************************************************

MAX_FIX  = 50
MAX_JUNC = 200
MAX_OBS  = 500


print "******************************************"
print
print "Least Square Leveling Network Adjustment"
print
print "******************************************"
print
print


/ V_fix_elv   = Vector of double to store FIX point elevation
/ V_fix_name  = Vector of String to store FIX point name
/ V_junc_name = Vector of String to store junction point name

V_fix_name = VecStr(MAX_FIX)
V_fix_elv = Vector(MAX_FIX)
V_junc_name = VecStr(MAX_JUNC)


fname     = String()
fname_out = String()
input "Input file name  : " fname
input "Output file name : " fname_out

Vstr_inp_tmp = VecStr()
Vstr_inp_tmp.load(fname)


n = Vstr_inp_tmp.size()

Vstr_inp = VecStr(n)

/ Re-read data again, trim left and right, and check no of fields
/ skip blank lines and comment lines


for i=0, n-1

   inp = Vstr_inp_tmp(i).ltrim().rtrim()

 / skip empty line or comment line
   if (inp.len()==0)
     continue
   end
   if (inp.left(1)=="/")
     continue
   end

   no = inp.no_word()

   if(no <>2 & no <>4)
     print "error in reading data : no of fields must be 2 or 4..."
     exit
   else
     Vstr_inp.pushback(inp)
   end
end




/ First read FIX points and elevations
/ Keep reading until a record that have 4 fields is found

n = Vstr_inp.size()

print
print "Reading FIX Points..."
set width 15

for i=0, n-1

   inp = Vstr_inp(i)

   no = inp.no_word()

   if(no==2)
     vec_tmp  = inp.extract(2)
     fix_name = vec_tmp(0)
     fix_elv  = vec_tmp(1).val()
     V_fix_name.pushback(fix_name)
     V_fix_elv.pushback(fix_elv)
     print fix_name;
     print fix_elv
   else
     break
   end
end

no_fix = V_fix_name.size()

/ Next read observation data
/ skip fix point data (having 2 fields)

print
print "Reading Junction Points..."
set width 15

no_obs = 0

for i=0, n-1

   inp = Vstr_inp(i)

   no = inp.no_word()

   if(no==2)
     continue
   else

      no_obs    = no_obs+1

      vec_tmp   = inp.extract(2)
      from_name = vec_tmp(0)
      to_name   = vec_tmp(1)

    / Check the from_name if it is a junction name
      ind1 = V_fix_name.find(from_name)
      if(ind1<0)
        ind11 = V_junc_name.find(from_name)
        if(ind11<0)
          print from_name
          V_junc_name.pushback(from_name)
        end
      end

    / Check the to_name if it is a junction name
      ind1 = V_fix_name.find(to_name)
      if(ind1<0)
        ind11 = V_junc_name.find(to_name)
        if(ind11<0)
          print to_name
          V_junc_name.pushback(to_name)
        end
      end

   end

end

no_junc = V_junc_name.size()

A  = Matrix(no_obs, no_junc)
L0 = Matrix(no_obs, 1)
Lb = Matrix(no_obs, 1)
P  = Matrix(no_obs, 1)



/ Now read observation data and create matrices A, Lb, L0

set width 15

no_obs = 0

for i=0, n-1

   inp = Vstr_inp(i)

   no = inp.no_word()

   if(no==2)
     continue
   else

     no_obs = no_obs+1

     vec_tmp   = inp.extract(4)
     from_name = vec_tmp(0)
     to_name   = vec_tmp(1)
     dist_km   = vec_tmp(2).val()
     dh        = vec_tmp(3).val()

     Lb(no_obs-1,0) = dh

   / P = weight and is inversely proportional to distance
     P (no_obs-1,0) = 1/dist_km

     ind1 = V_fix_name.find(from_name)
     if(ind1>=0)
       L0(no_obs-1,0) = -V_fix_elv(ind1)
     end
     ind1 = V_junc_name.find(from_name)
     if(ind1>=0)
       A(no_obs-1, ind1) = -1
     end

     ind1 = V_fix_name.find(to_name)
     if(ind1>=0)
       L0(no_obs-1,0) = V_fix_elv(ind1)
     end
     ind1 = V_junc_name.find(to_name)
     if(ind1>=0)
       A(no_obs-1, ind1) = 1
     end

   end

end

/ It's to do a computation
/ A = design matrix
/ Lb = Observation vector
/ La = Adjusted obs vector
/ X = Adjusted parameters
/ V = Residual

L = Lb - L0
N = A.tsp().multdiag(P)*A
y = A.tsp().multdiag(P)*L
X = N.inv()*y

V  = L - A*X
La  = Lb - V
var = V.tsp().multdiag(P)*V/(no_obs-no_junc)
sd  = sqrt(var(0,0))

/******************************************************************
/                     Writing Output file
/******************************************************************

fout_old = getfout()

set fout fname_out

fprint " ***********************************************"
fprint " * Least Square Leveling Network Adjustment *"
fprint " ***********************************************"
fprint
fprint
fprint
set width 4
fprint "Input file name   : " fname
fprint "Output file name  : " fname_out
fprint "no of observation : " int(no_obs)
fprint "no of parameter   : " no_junc
fprint "degree of freedom : " int(no_obs - no_junc)
fprint
fprint
fprint

fprint " List of Observation Data"
fprint
fprint " From To Dist dH "
fprint " Station Station (km) (m)"

fprint
for i=0, n-1

   inp = Vstr_inp(i)

   no = inp.no_word()

   if(no==2)
     continue
   else

     vec_tmp   = inp.extract(4)
     from_name = vec_tmp(0)
     to_name   = vec_tmp(1)
     dist_km   = vec_tmp(2).val()
     dh        = vec_tmp(3).val()

     set width 15
     set precision 4
     fprint from_name;
     fprint to_name;
     fprint dist_km;
     fprint dh

   end

end

fprint
fprint
fprint
fprint " List of Fix Elevation"
fprint
fprint " Name Elevation"
fprint

for i=0, no_fix-1
fprint V_fix_name(i);
fprint V_fix_elv(i)
end

fprint
fprint
fprint
fprint " Adjusted Data and Residual"
fprint
fprint " No Observed Adjusted Weight Residual Allowable error"
fprint " dH(m) dH(m) (1/km)  (m) (2nd order)"

fprint
for i=0, no_obs-1
   set width 5
   fprint int(i+1);
   set width 15
   fprint Lb(i,0);
   fprint La(i,0);
   fprint P(i,0);
   fprint V(i,0);
   fprint 0.0084*sqrt(1/P(i,0))
end

fprint
set precision 10
fprint " ***standard deviation of unit weight : " sd
fprint
fprint
fprint
fprint " Adjusted Elevation"
fprint
fprint " Name Elevation"
fprint
set precision 4

for i=0, no_junc-1
   print " ";
   fprint V_junc_name(i);
   fprint X(i,0)
end

fprint
fprint
fprint
fprint "Remark : - individual weight = 1/dist(KM)"
fprint " - individual sd = sqrt(1/weight)"
fprint " - allowable error(2nd order) = 0.0084 m * sqrt(dist KM)"
fprint " - Residual = Observed - Adjusted"
fprint
fprint " .............................."

/ restore fout
set fout fout_old
 

 

An example of input data file

Note :- data start with a list of FIX point and its elevation.  Each line must have 2 fields, namely station name and elevation.

         - It follows by a list of observation, height difference.  Each line must have 4 fields, namely from station name, to station name, distance (in KM), and observed height difference

         - Fields within a line are separated by at least one blank, space.  Station name must not have space in it.

         - Comment lines start with a sign "/", and empty lines are allowed in the file


/ List of FIX points
/ Name      Elevation

  BMP-402      126.0660
  BMP-404      129.0006

/ List of observation data
/ From           To     Dist(KM)    DH
  BMP-402     NEA_Y-06    6.712   -0.4968
  NEA_Y-06    lev-11      8.928   -0.9891
  lev-11      NEA_Y-15    3.249   -0.0292
  NEA_Y-15    lev-14      1.109   -0.6953
  lev-14      NEA_803     3.677    1.1463
  NEA_803     NEA_Y-28   11.083    0.1726
  NEA_Y-28    BMP-404    48.808    3.8404
  lev-11      NEA_Y-15   15.446   -0.0400
  NEA_Y-28    NEA_Y-01   10.866    0.3917
  NEA_Y-01    lev-34      3.910    0.7717
  lev-34      NEA_Y-06    7.190   -0.7839
  lev-34      NEA_Y-20    3.663    0.2539
  NEA_Y-20    NEA_Y-01    7.292   -1.0506
  NEA_Y-20    lev-14      6.585   -2.7479
  BMP-404     lev-134     5.708   -2.7177
  lev-134     NEA_Y-48   16.556    0.8032
  NEA_Y-48    NEA_Y-33   14.493   -1.7467
  NEA_Y-33    AZ_S.13    12.588    2.7237
  AZ_S.13     NEA_803     5.816   -3.0368
  NEA_Y-33    lev-78      3.658    1.8177
  lev-78      AZ_S.13    15.241    0.8667
  lev-78      lev-106    24.920    1.2312
  lev-106     NEA_Y-48    1.485   -1.2787
  lev-106    lev-134 3    0.123   -2.1132

 

An example of output file

                            
      ***********************************************
      * Least Square Levelling Network Adjustment *
      ***********************************************



Input file name   : ytlslev.txt
Output file name  : ytlslev_out.txt
no of observation : 24
no of parameter   : 15
degree of freedom : 9



                List of Observation Data

  From            To                    Dist           dH
Station         Station                 (km)           (m)

BMP-402        NEA_Y-06                6.7120        -0.4968
NEA_Y-06       lev-11                  8.9280        -0.9891
lev-11         NEA_Y-15                3.2490        -0.0292
NEA_Y-15       lev-14                  1.1090        -0.6953
lev-14         NEA_803                 3.6770         1.1463
NEA_803        NEA_Y-28               11.0830         0.1726
NEA_Y-28       BMP-404                48.8080         3.8404
lev-11         NEA_Y-15               15.4460        -0.0400
NEA_Y-28       NEA_Y-01               10.8660         0.3917
NEA_Y-01       lev-34                  3.9100         0.7717
lev-34         NEA_Y-06                7.1900        -0.7839
lev-34         NEA_Y-20                3.6630         0.2539
NEA_Y-20       NEA_Y-01                7.2920        -1.0506
NEA_Y-20       lev-14                  6.5850        -2.7479
BMP-404        lev-134                 5.7080        -2.7177
lev-134        NEA_Y-48               16.5560         0.8032
NEA_Y-48       NEA_Y-33               14.4930        -1.7467
NEA_Y-33       AZ_S.13                12.5880         2.7237
AZ_S.13        NEA_803                 5.8160        -3.0368
NEA_Y-33       lev-78                  3.6580         1.8177
lev-78         AZ_S.13                15.2410         0.8667
lev-78         lev-106                24.9200         1.2312
lev-106        NEA_Y-48                1.4850        -1.2787
lev-106        lev-134                30.1230        -2.1132



         List of Fix Elevation

  Name                 Elevation

BMP-402                126.0660
BMP-404                129.0006



                     Adjusted Data and Residual

  No      Observed      Adjusted        Weight      Residual  Allowable error
            dH(m)         dH(m)        (1/km)          (m)      (2nd order)

   1       -0.4968       -0.4977        0.1490        0.0009       0.0218
   2       -0.9891       -0.9857        0.1120       -0.0034       0.0251
   3       -0.0292       -0.0301        0.3078        0.0009       0.0151
   4       -0.6953       -0.6949        0.9017       -0.0004       0.0088
   5        1.1463        1.1466        0.2720       -0.0003       0.0161
   6        0.1726        0.1713        0.0902        0.0013       0.0280
   7        3.8404        3.8251        0.0205        0.0153       0.0587
   8       -0.0400       -0.0301        0.0647       -0.0099       0.0330
   9        0.3917        0.3939        0.0920       -0.0022       0.0277
  10        0.7717        0.7791        0.2558       -0.0074       0.0166
  11       -0.7839       -0.7802        0.1391       -0.0037       0.0225
  12        0.2539        0.2590        0.2730       -0.0051      0.0161
  13       -1.0506       -1.0382        0.1371       -0.0124       0.0227
  14       -2.7479       -2.7499        0.1519        0.0020       0.0216
  15       -2.7177       -2.7187        0.1752        0.0010       0.0201
  16        0.8032        0.8117        0.0604       -0.0085       0.0342
  17       -1.7467       -1.7573        0.0690        0.0106       0.0320
  18        2.7237        2.7058        0.0794        0.0179       0.0298
  19       -3.0368       -3.0379        0.1719        0.0011       0.0203
  20        1.8177        1.8202        0.2734       -0.0025       0.0161
  21        0.8667        0.8856        0.0656       -0.0189       0.0328
  22        1.2312        1.2176        0.0401        0.0136       0.0419
  23       -1.2787       -1.2805        0.6734        0.0018       0.0102
  24       -2.1132       -2.0922        0.0332       -0.0210       0.0461

   ***standrad deviation of unit weight : 0.0040270687



     Adjusted Elevation

    Name        Elevation

  NEA_Y-06       125.5683
  lev-11         124.5826
  NEA_Y-15       124.5525
  lev-14         123.8576
  NEA_803        125.0042
  NEA_Y-28       125.1755
  NEA_Y-01       125.5694
  lev-34         126.3486
  NEA_Y-20       126.6076
  lev-134        126.2819
  NEA_Y-48       127.0935
  NEA_Y-33       125.3363
  AZ_S.13        128.0421
  lev-78         127.1565
  lev-106        128.3741



Remark : - individual weight          = 1/dist(KM)
         - individual sd              = sqrt(1/weight)
         - allowable error(2nd order) = 0.0084 m * sqrt(dist KM)
         - Residual                   = Observed - Adjusted

              ..............................

 

| HOME |   BACK   |

Mastermind.prg

 

/****************************************************************************
/                              MASTERMIND GAME
/****************************************************************************

/ This program is to play the "MASTERMIND" game.
/ Instruction:
/      on running the program, the user is asked
/      to guess a series of FOUR SECRET NUMBERS
/
/ by doing this, the user must key in 4 numbers (0-9)
/ separated by at least one space.

/ Then the program will compare your numbers with the ones stored in the program
/ and report that how many numbers are in the correct places, and
/ how many number that are correct but in the wrong position.

/ Author : goodgirl
/ Date   : 7 March 2002

set randmax 5000

MAX_NO = 9
NO     = 4

/ vec_code is to store the code numbers, secret numbers
vec_code = VecInt(4)
for i = 0, NO-1
   a = rand()%(MAX_NO+1)
   vec_code.pushback(a)
end

print
print
print " ********************"
print "    Mastermind Game"
print " *******************"
print
print
print "Keep on inputing 4 numbers (0-9) separated by a space"
print "To quit the program, press CTRL C followed by <Enter>"
print
print

temp = 1
no_guess = 0

while temp== 1
   no_guess = no_guess + 1
   s = string()
   input "the four numbers : " s
   n = s.no_word()
   if n<>NO
      print "error : must input exactly 4 numbers, separated by a space"
      continue
   end

   vec_guess = s.extract(NO)
   vec_guess_tmp = VecInt(NO)

   / put the guess numbers in a vector
   for i=0, NO-1
      a = round(vec_guess(i).val())
      vec_guess_tmp.pushback(a)
   end

   / no_cor_pos = correct number and correct position
   / no_cor_num = correct number but wrong position

   no_cor_pos = 0
   no_cor_num = 0

   vec_code_tmp = vec_code

   wrong_number = 99
   / First check if there are correct number + position
   for i=0, NO-1
      a = vec_guess_tmp(i)
      b = vec_code_tmp(i)
      / If found, alter the code number and the guess number
      / so that they cannot involve in later processing
      / also keep increasing value of wrong_number,
      / so that it is always unique
      if (a==b)
         no_cor_pos       = no_cor_pos + 1
         vec_code_tmp(i)  = wrong_number
         vec_guess_tmp(i) = wrong_number + 1
         wrong_number     = wrong_number + 2
      end
   end

   / Now check the guess number one by one, against the code numbers
   for j = 0,NO-1
      a = vec_guess_tmp(j)
      for k = 0, NO-1
         b = vec_code_tmp(k)
         if (a==b)
            no_cor_num      = no_cor_num + 1
            vec_code_tmp(k) = wrong_number

            wrong_number    = wrong_number + 1
            break
         end
      end
   end

   set width 3

   print " No : " int(no_guess);
   print " correct position : " int(no_cor_pos);
   print " correct number (wrong position) : " int(no_cor_num)

   if (no_cor_pos==4)
     print
     print "*** Congratulation!!! the numbers are : " s
     exit
   end
end


 

 

 


| HOME |   BACK   |