| HOME |
BACK |
LDT.fun
/****************************************************************************
/
LDT
/****************************************************************************
/ This function is to calculate Lens distortion at an x y photo coordinate
/ of a Photo object
/ Author : Parallax
/ Date : June,
2017
Function Pt2D lens_distortion
= ldt(Photo a, double x, double y)
no_para = a.io_para12().size()
if(no_para <> 4 & no_para <> 6 & no_para <> 7 & no_para <> 12)
print "Error IO parameters(forward) must be 4 or 6 or 7 or 12 values"
return
end
// analog photo
if (no_para ==4 | no_para==6 |
no_para==12)
lens_distortion = -999
return
// digital
photo
else
rdist2 = x^2 + y^2
dx_rd =
a.io_para12()(0)*rdist2*x + a.io_para12()(3)*rdist2*rdist2*x
dy_rd =
a.io_para12()(0)*rdist2*y + a.io_para12()(3)*rdist2*rdist2*y
dx_dd
= a.io_para12()(2)*(rdist2+2.0*x*x) + 2.0*a.io_para12()(5)*x*y
dy_dd =
2.0*a.io_para12()(2)*x*y + a.io_para12()(5)*(rdist2-2.0*y*y)
dx_af
= -a.io_para12()(1)*x + a.io_para12()(4)*y
dy_af = a.io_para12()(4)*x
end
lens_distortion = Pt2D(dx_rd + dx_dd + dx_af, dy_rd +
dy_dd + dy_af)
return
|
| 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 |
|
