(ฟรี)โปรแกรมภาษาไพทอนบนเครื่องคิดเลขคาสิโอ fx-9750GIII fx-9860GIII และ fx-cg50 โปรแกรมพื้นฐานงานสำรวจชุดที่ 4 (COGO Selected Serie 4)

ตอนนี้มาถึงโปรแกรมพื้นฐานงานสำรวจชุดที่ 4 (COGO Selected Serie 4) สำหรับเครื่องคิดเลขคาสิโอ fx-9750GIII, fx-9860GIII และ fx-cg50 PRIZM สามรุ่นที่รองรับภาษาไพทอนหรือไมโครไพทอน ได้ในขณะนี้ หาซื้อได้ในเมืองไทย ราคาย่อมเยาที่สุดคือ fx-9750GIII ที่ราคาประมาณสามพันบาท ถ้ามีงานการทำเป็นหลักเป็นแหล่งแล้วไม่น่าแพง บางทีเราซื้อโทรศัพท์มือถือได้ราคาเป็นเรือนหมื่นไม่คิดอะไรมาก แต่กับเครื่องคิดเลขคิดแล้วคิดอีก สำหรับผมแล้วราคาขนาดนี้ถือว่าคุ้มค่ามาก ผมไปนั่งพลิกดูด้านหลังเครื่อง fx-9860GIII ปรากฎเห็น “Made in Thailand” ค่อนข้างประหลาดใจพอสมควร พอไปพลิกดูเครื่อง fx-cg50 กลับเป็น “Made in China”

เครื่องคิดเลขคาสิโอ fx-9750GIII, fx-9860GIII และ fx-CG50 PRIZM

วิธีการติดตั้งโปรแกรมไพทอนลงบนเครื่องคิดเลข

ผมเคยเห็นแผ่นผีลงโปรแกรมคำนวณวงรอบ “Traverse Pro” ที่ศูนย์ไอทีหลายแห่งเมื่อประมาณ 7-8 ปีที่แล้ว ก็สงสัยอยู่ครามครันว่าโปรแกรมมันแจกฟรีในบล็อกผมอยู่แล้ว ทำไมคนยังต้องซื้อแผ่นผีอีก ก็มานึกได้ว่าในยุคนั้นสังคมโซเชียลยังไม่เบ่งบานเท่าปัจจุบัน คนที่ไม่รู้ก็ไม่มีเครือข่ายบอกกันต่อ เช่นเดียวกัน ผมคิดว่าการก็อปปี้โปรแกรมภาษาไพทอนลงเครื่องคิดเลขเป็นเรื่องง่ายๆ ที่ใครๆเข้าใจได้ คงจะไม่ใช่ เลยจัดทำวิธีการติดตั้งลงโปรแกรมภาษาไพทอนลงเครื่องคิดเลข ขั้นตอนก็ไปดาวน์โหลดโปรแกรมจากบล็อกนี้เข้าเครื่องคอมพิวเตอร์ จากนั้นเอาสายเคเบิ้ลมินิยูเอสบี มาเสียบเครื่องคอมพิวเตอร์และเครื่องคิดเลขจากนั้นทำตามคู่มือลิ๊งค์นี้ครับ

ชุดโปรแกรมพื้นฐานงานสำรวจ COGO (Coordinate Geometry) ชุดที่ 4

มาถึงโปรแกรมพื้นฐานงานสำรวจชุดที่ 4 (COGO Selected Serie 4) ซี่รี่ย์ว่าด้วยเรื่องการยืดจุดเส้นตรงในสามมิติไปยังระนาบ (plane) หรือการยืดจุดจากเส้นตรงสามมิติเมื่อทราบระยะทาง ตลอดจนการคำนวณหาพื้นที่ (area) และสุดท้ายคือโปรแกรมมาร์คแกน (mark axis) ซึ่งรายละเอียดมาดูกันอีกทีว่าใช้งานและการประยุกต์ใช้งานอย่างไรบ้าง

โปรแกรมพื้นฐานงานสำรวจชุดที่ 4 บนเครื่องคิดเลขคาสิโอ fx-9750GIII

ส่วนประกอบของโปรแกรม

สำหรับโปรแกรมพื้นฐานงานสำรวจในชุดนี้จะจัดโปรแกรมย่อยเล็กๆ ไว้ 4 โปรแกรม เนื่องจากหน้ากว้างของหน้าจอเครื่องคิดเลขมีแค่ 21 ตัวอักษร ดังนั้นบางอย่างผมใช้คำย่อเช่น Ext ย่อมาจาก Extend แปลว่ายืดหรือขยาย และคำว่า Pt ย่อมาจาก point คือจุดนั่นเอง

  1. Ext Line to Plane เมื่อรู้ค่าพิกัดของเส้นตรง (X,Y,Z) 2 จุด และระนาบที่กำหนดด้วยจุด (X,Y,Z) 3 จุด สามารถยืดเส้นตรงไปแตะระนาบแล้วหาจุดพิกัด (X,Y,Z) นั้นได้
  2. Ext Pt from 3D Line เมื่อรู้ค่าพิกัดของเส้นตรง (X,Y,Z) 2 จุด ต้องการยืดระยะทางไปตามที่กำหนด สามารถหาจุดพิกัด (X,Y,Z) นั้นได้
  3. Calculate Area เมื่อกำหนดจุดค่าพิกัดของแปลงพื้นที่ดินรอบรูป สามารถป้อนค่าพิกัดเพื่อคำนวณหาพื้นที่ได้
  4. Mark Axis ต้องการมาร์คแกนโดยที่รู้ค่า design และที่หน้างานทำการวางจุดด้วยกล้องโททัลสเตชั่น ทำการคำนวณด้วยเครื่องคิดเลขเพื่อหาจุดได้ 2 จุดลากเส้นเชื่อม เมื่อครบ 4 จุด จะได้เส้นตรงสองเส้นตัดกันเป็นมุม 90 องศา

วิธีการใช้งานโปรแกรม

กดคีย์ “MENU” ที่เครื่องคิดเลข กดคีย์ลูกศรเลื่อนไปที่ไอคอน “Python” กดเข้าไป ตรงนี้ผมสร้างโฟลเดอร์ชื่อ “PROGRAM” ผมเอาซอร์สโค้ดโปรแกรมไพทอนมาไว้ที่นี่ เมื่อกด “EXE” เข้าไปจะเห็นไฟล์ซอร์สโค้ดไพทอนเรียงกันประมาณนี้ ตอนนี้จะรันไฟล์ COGOSSE4.py ซึ่งจะต้องมีไฟล์ไลบรารี pbrutils.py ด้วย

Ext Line to Plane

โปรแกรมช่วยคำนวณการยืดเส้นตรงในระบบสามมิติไปตัดหรือทะลุระนาบสามมิติ โดยที่ทราบค่าพิกัดของเส้นตรง 2 จุด (X, Y, Z) และทราบค่าพิกัดของระนาบ 3 จุด (X, Y, Z)

กำหนดค่าพิกัดเส้นตรงจุดที่ 1 (N=5002.002, E=3001.001, Elev=49.090) และจุดที่ 2 (N=4996.288, E=3005.543mE, Elev=54.190) และกำหนดระนาบ

  1. จุดที่ 1 N=5010.456, E=3010.123, Elev=55.678 m
  2. จุดที่ 2 N=5008.939, E=3005.515, Elev=52.010 m
  3. จุดที่ 3 N=5009.999, E=3007.777, Elev=58.888 m

ต้องการหาจุดที่เส้นตรงไปตัดระนาบและหามุมที่ทำกับระนาบ

ที่เมนูกดคีย์ “1” และกดคีย์ “EXE” โปรแกรมจะถามค่าพิกัดเส้นตรงจุดที่ 1 และจุดที่ 2 จากนั้นจะถามค่าพิกัดจุดของระนาบจำนวน 3 จุด

จะได้จุดตัด N = 5006.301, E = 2997.584 และ Z = 45.253 และมุมที่เส้นตรงทำกับ plane 52.0219724 องศาจากระนาบ X-Y

Ext Pt from 3D Line

โปรแกรมสำหรับยืดเส้นตรงในระบบ 3 มิติไปตามระยะทางที่กำหนด จากไดอะแกรมด้านล่างรูปซ้ายเป็นเป้าที่มี 2 ปริซึม กล้องโททัลสเตชัน สามารถวัดค่าพิกัดจุด 1 (N1, E1, Z1) และจุดที่ 2 (N2, E2, Z2) เมื่อทราบความสูงโพล d ก็สามารถหาค่าพิกัดของจุด P ที่ซ่อนอยู่ได้

ส่วนรูปด้านขวาจะเป็นจุด P ที่มุมห้องที่กล้องโททัลสเตชันมองไม่เห็น ลากเทปเหล็กออกมาวัดค่าพิกัดแบบ Reflector-less ที่ 4 เมตร และ 2 เมตร ดังนั้นระยะทางที่ยืดไปหามุมห้องเท่ากับ 2 เมตร สามารถหาค่าพิกัดที่มุมห้องได้เช่นกัน

โจทย์กำหนดเส้นตรง จุดที่ 1 (N = 5184.542, E = 1723.884, Elev = 68.888) และจุดที่ 2 (N = 5185.409, E = 1732.676, Elev = 69.010) ยืดจุดจากจุดที่ 2 ไปอีก 3.075 เมตร ต้องการทราบค่าพิกัดของจุดนี้

ที่เมนูกดคีย์ “2” และกดคีย์ “EXE” โปรแกรมจะถามค่าพิกัดเส้นตรงจุดที่ 1 และจุดที่ 2 จากนั้นจะถามระยะทางที่ต้องการยืดหรือขยายออกไป

จะได้คำตอบคือ N = 5185.711, E = 1735.736 และ Elev = 69.052

Calculate Area

การคำนวณพื้นที่เมื่อกำหนดค่าพิกัดฉากของจุดรอบแปลงปิด โดยการใช้สูตร Shoelace formula หรือสูตรผูกเชือกรองเท้า ที่อาศัยการคูณไขว้ ได้เท่าไหร่จับมาลบกันแล้วหารด้วยสอง

ไดอะแกรมแสดงการคำนวณพื้นที่โดยใช้สูตรคูนไขว้หรือสูตร shoelace formula

กำหนดจุดตารางด้านล่าง ต้องการหาพื้นที่

PointNorthingEasting
1565.323m544.145m
2564.690m577.134m
3587.898m579.579m
4594.486m565.478m
5585.709m533.921m

ที่เมนูกดคีย์ “3” และกดคีย์ “EXE” โปรแกรมจะถามค่าพิกัดจุดที่ 1 จุดที่ 2 จุดที่ 3 เมื่อต้องการปิดให้ย้อนมาป้อนจุดที่ 1 อีกครั้ง (รูปด้านล่างเน้นด้วยสี่เหลี่ยมสีแดง) โปรแกรมจะหยุดให้ป้อนค่าพิกัดและคำนวณพื้นที่มาให้

ได้คำตอบคือพื้นที่ 1025.991 ตารางเมตร

Mark Axis

ในงานก่อสร้างบางอย่างเช่นหัวเสาหรือฟุตติ้ง จำเป็นต้องมีการตีแกนให้โดยการตีเต๊า เพื่อเป็นเส้น reference ให้ทางช่างก่อสร้างหน้างานสามารถวัดออฟเซ็ทด้วยตลับเมตรหรือเรียกง่ายๆว่าแทงตลับเมตร เพื่อติดตั้งฟอร์มเวิร์คหรือทำอย่างอื่นได้ง่ายๆ สิ่งที่ช่างสำรวจต้องรู้คือค่าพิกัดจุดศูนย์กลางและอะซิมัทของแกน รูปด้านล่างเป็นหัวเสากว้าง 1 เมตร ยาว 2 เมตร จุด C หรือจุดศูนย์กลางอยู่ที่ N = 2000m, E= 1000m อะซิมัทของหัวเสาเอียงอยู่ 40°12’30” จากทิศเหนือ โจทย์ต้องการมาร์คจุดศูนย์กลาง และมาร์คจุดหมายเลข 1, 2, 3 และ 4 จากนั้นก็ตีเต๊าเส้น 1-3 (แกน Y) เส้น 2-4 (แกน X)

ถ้าคิดว่าแกน YX มีจุดกำเนิด (origin) ผ่านจุดศูนย์กลางหัวเสาพอดี ค่าพิกัดจุดที่ 1 จะเป็น Y = 2, X = 0 จุดที่ 2 Y = 0, X = 1 จุดที่ 3 Y = -2, X = 0 และจุดที่ 4 Y = 0, X = -1

ที่เมนูกดคีย์ “4” และกดคีย์ “EXE” โปรแกรมจะถามค่าพิกัด (Design) จากนั้นจะถามค่าพิกัดที่เก็บได้ (Actual) จากหน้างาน ณ ตอนนั้นจากกล้องโททัลสเตชั่น จากนั้นจะคำนวณค่า Deviated เมื่อเทียบกับแกน NE และ แกน YX เราจะใช้ค่านี้ในการขยับจุดเพื่อให้เข้าหาแกน

ต้องการมาร์คจุด C หรือจุดศูนย์กลาง วัดค่าพิกัดได้ N = 1999.982 E = 1000.025

จะได้ผลลัพธ์ NE: -0.018 0.025 และ YX: 0.002 0.031 ตัวเลขนี้บอกอะไรเราได้บ้าง ผมขยายรูปด้านบนและเขียนขนาดมิติ ค่า ΔN = -0.018 ΔE = +0.025 เมื่อเทียบกับพิกัดฉากที่เราใช้คือ NE ส่วนพิกัดฉาก YX คือพิกัดฉากของหัวเสาโดยที่ค่า ΔY = +0.002 ΔX = +0.031แสดงว่าจุดใหม่ที่จะวางเกินไป 0.002 ต้องขยับตามแกน Y ลงมาให้เท่ากับศูนย์ และเช่นเดียวกันเกินไป 0.031 ต้องขยับตามแกน X ไปด้านซ้ายให้เท่ากับศูนย์ ขยับจุดใหม่ป้อนค่าพิกัดเข้าโปรแกรมปรับจนได้ค่าเบี่ยงเบนเท่ากับศูนย์ YX: 0.000 0.000 จะได้จุดศูนย์กลาง

ต่อไปมาร์คจุดที่ 1 ถ้าคิดว่าทำการวางจุดอ่านจากกล้องโททัลสเตชั่นได้ N = 2001.514 E = 1001.292 ทำการป้อนค่าเข้าไปในโปรแกรม

จะได้ค่าเบี่ยงเบน NE: 1.514 1.292 และ YX: 1.990 0.009 ในที่นี้เราต้องการค่าพิกัดเทียบกับแกน YX โดยที่ Y = 2, X = 0 (ดูรูปแรก) ดังนั้นต้องขยับไปตามแกน Y ขี้นไปเท่ากับ 0.010 ถึงจะได้ 2 เมตรพอดี และตัวเลข 0.009 เกินจากแกน X ไปด้านขวา ต้องขยับไปด้านซ้าย 0.009 เพื่อลดลงเท่ากับศูนย์ เมื่อได้จุดที่ต้องการจะได้ค่า YX: 2.000 0.000 ทำการวางจุดที่ 2, 3 และ 4 ด้วยวิธีการเดียวกันนี้ จากนั้นทำการตีเต๊าเชื่อมจุดจะได้เส้น 2 เส้นตัดกันเป็นมุมฉาก

สรุป

ก็จบไปแล้วครับชุดโปรแกรมพื้นฐานงานสำรวจชุดที่ 4 (COGO Selected Serie 4) จะเห็นแต่ละโปรแกรมไม่มีกราฟฟิควาดรูปมาให้ดูด้วย ผมพยายามลองวาดด้วยไลบรารี matplotl.py แต่โปรแกรมหลุดบ่อย เลยล้มเลิกไป อนาคตถ้าคาสิโอเตรียมไลบรารีมาดีและเสถียรกว่านี้ ผมจะย้อนมาแก้ไขโค้ด เพิ่มส่วนวาดรูป

หลายปีที่ผ่านมายอมรับว่าไม่ได้ติดตามตลาดเครื่องคิดเลขมากนัก ผมลองดูเครื่องคิดเลขในท้องตลาดตอนนี้ยี่ห้อ TI ก็มีหลายรุ่นที่เขียนภาษาไพทอนได้ แต่ที่เหนือกว่าคาสิโอ คือมีไลบรารีต่างๆที่ TI เตรียมมาให้ด้วยเช่นฟังก์ชั่นการวาดรูป การอ่านเขียนไฟล์ แต่เมื่อข้ามไปยี่ห้อ HP รุ่น Prime พบว่าสถาปัตยกรรมเหมือนของโทรศัพท์มือถือคือใช้ซีพียูของอาร์ม (Arm V7) ซีพียูแรงกว่าแต่ต้องแลกกับพลังงานที่มากกว่า ใช้ถ่านชาร์จอาจจะต้องชาร์จกันบ่อย ที่น่าสนใจคือโปรแกรมมิ่งใช้ภาษาปาสคาล (Pascal-like) ที่ผมคุ้นเคยอยู่ ว่างๆว่าจะถอยมาเล่นๆอีกสักเครื่อง ราคาก็สูงอยู่ขนหน้าแข้งคงจะร่วงอีกหลายเส้น โปรดติดตามบทความตอนต่อไปครับ

ซอร์สโค้ดโปรแกรมพื้นฐานงานสำรวจชุดที่ 4 (COGOSSE4.py)

from pbrutils import split_angle,deg2dms,DMS2str
from pbrutils import calcAzimuth,calcDistance,calcDiff2Angle
from pbrutils import calcNextAzimuth
from pbrutils import calcIntersectionAzimuthAzimuth,equation_plane
from pbrutils import PI,DEG2RAD,RAD2DEG,mean,calcNextCoordinate2D
from pbrutils import isect_line_plane_v3_4d
from pbrutils import calcExtend3DLineCoords,area_by_shoelace
import math

def print_menu():
  print('COGO Selected Serie 4')
  print('1:Ext Line to Plane')
  print('2:Ext Pt from 3D Line')
  print('3:Calculate Area')
  print('4:Mark Axis')
  print('5:Exit')

def extend_line_to_plane():
  print('Extend Line to Plane:')
  print("Input Line:")
  #n1,e1,z1=1.0,4.0,1.0
  n1=float(input("N1=?"))
  e1=float(input("E1=?"))
  z1=float(input("Z1=?"))
  #n2,e2,z2=-2.0,8.0,-2.0
  n2=float(input("N2=?"))
  e2=float(input("E2=?"))
  z2=float(input("Z2=?"))
  print("Input 3 pts of Plane:")
  #n3,e3,z3=2.0,-1.0,4.0
  n3=float(input("N3=?"))
  e3=float(input("E3=?"))
  z3=float(input("Z3=?"))
  #n4,e4,z4=1.0,2.0,3.0
  n4=float(input("N4=?"))
  e4=float(input("E4=?"))
  z4=float(input("Z4=?"))
  #n5,e5,z5=3.0,1.0,2.0
  n5=float(input("N5=?"))
  e5=float(input("E5=?"))
  z5=float(input("Z5=?"))
  a,b,c,d=equation_plane(n3,e3,z3,n4,e4,z4,n5,e5,z5)
  plane=[a,b,c,d]
  pi,ang=isect_line_plane_v3_4d([e1,n1,z1],[e2,n2,z2],plane)
  if pi:
    print("Intersection point:")
    print("N= {0:.3f}".format(pi[1]))
    print("E= {0:.3f}".format(pi[0]))
    print("Z= {0:.3f}".format(pi[2]))
    print("Angle: {:.7f} deg".format(ang*RAD2DEG))
  else:
    print("Not intersect plane!")
  input("Press any key")  

def extend_point_from_3Dline():
  print("Ext point from 3D line.")
  print("Input Line:")
  n1=float(input("N1=?"))
  e1=float(input("E1=?"))
  z1=float(input("Z1-?"))
  n2=float(input("N2=?"))
  e2=float(input("E2=?"))
  z2=float(input("Z2=?"))
  print("Input extend distance:")
  ext=float(input("Extend dist=?"))
  n3,e3,z3=calcExtend3DLineCoords(n1,e1,z1,n2,e2,z2,ext)
  print("Ext pt coordinates:")
  print("N= {:.3f}".format(n3))
  print("E= {:.3f}".format(e3))
  print("Z= {:.3f}".format(z3))
  input("Press any key") 

def calculate_area():
  pts=[]
  i=1
  print("Calculate area.")
  print("Input points:")
  n0=float(input("N1=?" ))
  e0=float(input("E1=?" ))
  pts.append((n0,e0))
  x0,y0=e0,n0
  f = True
  while f:
    y=float(input("N{}=?".format(i+1)))
    x=float(input("E{}=?".format(i+1)))
    pts.append((y,x))
    if (i>2):
      if (x==x0) and (y==y0):
        f=False
    i+=1    
  y,x=zip(*pts)   
  area=area_by_shoelace(x,y)
  print("Result of calculation:")
  print("Area= {:0.3f}".format(area))
  input("Press any key") 

def mark_axis():
  f=True
  print("Mark Axis.")
  print("Input Design:")
  nde=float(input("N=?" ))
  ede=float(input("E=?" ))
  az1=input("Azimuth=?")
  azi1=split_angle(az1)*DEG2RAD 
  while f:
    print("Input Actual:")
    nac=float(input("N=?" ))
    eac=float(input("E=?" )) 
    azide2ac=calcAzimuth(nde,ede,nac,eac)
    dist=calcDistance(nde,ede,nac,eac)
    print("Dist= {:.3f}".format(dist))
    print("Deviated:")
    print("NE: {:.3f} {:.3f}".format(nac-nde,eac-ede))
    xy1=dist*math.cos(azide2ac-azi1)
    xy2=dist*math.sin(azide2ac-azi1)
    print("YX: {:.3f} {:.3f}".format(xy1,xy2))
    print("Total {:.3f}".format(dist))
    con=input("Continue (Y/N)?")
    if con=="Y" or con=="y":
      f=True
    else:
      f=False

loop=True
choice=-1
ex_choice=-1
while loop:      
  print_menu()
  choice=int(input('Selection[1-5]'))  
  if (choice==5):
    loop=False
  elif (choice==1):
    loop=True
    extend_line_to_plane()
    ex_choice=1
  elif (choice==2):
    loop=True
    extend_point_from_3Dline()
    ex_choice=2
  elif (choice==3):
    calculate_area()
    loop=True
    ex_choice=3
  elif (choice==4):
    loop=True
    mark_axis()
    ex_choice=4

ไลบรารี pbrutils.py

import math
PI=math.pi
DEG2RAD=PI/180.0
RAD2DEG=180.0/PI
TOLERANCE=1e-10
def clear():
  for i in range(6):
    print()
	
def split_angle(st):
  ss=st.split('-')
  j=0
  d,m,s=0.0,0.0,0.0
  for c in ss:
    if j==0:
      d=float(c)
    elif j==1:
      m=float(c)
    elif j==2:
      s=float(c)
    j+=1    
  return d+m/60.0+s/3600.0

def deg2dms(dd):
  sign=1
  if (dd<0):
    sign=-1
  dd=abs(dd)
  minutes,seconds=divmod(dd*3600,60)
  degrees,minutes=divmod(minutes,60)
  return (sign,degrees,minutes,seconds)

def DMS2str(degree,minute,second,numdec):
  degree=abs(degree); minute=abs(minute); second=abs(second)
  s='{:.%df}' % numdec
  ss=s.format(second)
  smin=s.format(60.0)
  mm='{:.0f}'.format(minute)
  if (ss==smin):
    minute+=1
    ss=s.format(0.0)
    if (minute>=60):
      mm='0'
      degree+=1
    else:
      mm='{:.0f}'.format(minute)
  return '{:.0f}'.format(degree)+"-"+mm+"-"+ss

def calcDistance(n1,e1,n2,e2):
  return (math.sqrt(math.pow(n2-n1, 2.0) + math.pow(e2-e1, 2)))

def calcAzimuth(n1,e1,n2,e2):
  """ atan2 not same as excel function"""
  angle=math.atan2(e2-e1,n2-n1)
  if angle>=0:
    return angle
  else:
    return angle+2*PI
    
def calcDiff2Angle(azi2,azi1):
  temp=azi2-azi1
  if (temp<0):
    return temp + 2*PI
  else:
    return temp
	
def calcNextAzimuth(currentang,prevazi):
  temp=currentang+prevazi
  if (temp<PI):
    temp=temp+PI
  elif (temp>3*PI):
    temp=temp-3*PI
  else:
    temp=temp-PI
  return temp

def calcNextCoordinate2D(n,e,prevazi,horang,griddist):
  nextazi=calcNextAzimuth(horang,prevazi)
  nt=n+math.cos(nextazi)*griddist
  et=e+math.sin(nextazi)*griddist
  return [nt,et]

def calcNextCoordinate3D(n,e,z,prevazi,horang,verang,slopedist,hi,ht):
  hordist=slopedist*math.sin(verang)
  verdist=slopedist*math.cos(verang)
  nextazi=calcNextAzimuth(horang,prevazi)
  nt=n+math.cos(nextazi)*hordist
  et=e+math.sin(nextazi)*hordist
  zt=z+hi+verdist-ht
  return [nt,et,zt]

def mean(num):
  sum_num = 0
  for t in num:
    sum_num = sum_num + t           
  avg = sum_num / len(num)
  return avg

def calcIntersectionFrom4Points(pt1, pt2, pt3, pt4):
    '''Compute the intersection points which given 4 points.
       If intersection on line return point otherwise return None.
    '''
    azimuth12 = calcAzimuth(pt1, pt2)
    distance12 = calcDistance(pt1, pt2)
    azimuth34 = calcAzimuth(pt3, pt4)
    distance34 = calcDistance(pt3, pt4)
    #print('++++ pt1 pt2 pt3 pt4 azimuth12 azimuth34 = ', pt1, pt2, pt3, pt4, azimuth12, azimuth34)

    if (math.abs(azimuth12 - (0.5 * PI)) <= TOLERANCE) or (math.abs(azimuth12 - (1.5 * PI)) <= TOLERANCE):
        pti = [(pt1[1] - pt3[1]) * math.tan(azimuth34) + pt3[0], pt1[1]]
    elif (azimuth34 == (0.5 * PI)) or (azimuth34 == (1.5 * PI)):
        pti = [(pt3[1] - pt1[1]) * math.tan(azimuth12) + pt1[0], pt3[1]]
    elif (azimuth34 == azimuth12):
        return []
    else:
        y = (math.tan(azimuth34) * pt3[1] - math.tan(azimuth12) * pt1[1] + pt1[0] 
            - pt3[0]) / (math.tan(azimuth34) - math.tan(azimuth12))
        x = (y - pt1[1]) * math.tan(azimuth12) + pt1[0]
        pti = [x, y]
    
    dist1ni = calcDistance(pt1, pti)
    dist2ni = calcDistance(pt2, pti)
    dist3ni = calcDistance(pt3, pti)
    dist4ni = calcDistance(pt4, pti)
    if (math.abs((dist1ni + dist2ni) - distance12) <= TOLERANCE) and (math.abs((dist3ni + dist4ni) - distance34) <= TOLERANCE):
        return pti
    else:
        return []
    
def calcIntersectionAzimuthAzimuth(pt1,azimuth1,pt2,azimuth2):
  '''Azimuth is radian'''
  n1=pt1[0] #N
  e1=pt1[1] #E
  n2=pt2[0]
  e2=pt2[1]
  delty=n2-n1
  deltx=e2-e1
  ni=(math.tan(azimuth2)*n2-math.tan(azimuth1)*n1+e1-e2)/(math.tan(azimuth2)-
     math.tan(azimuth1))
  ei=(ni-n1)*math.tan(azimuth1)+e1
  return [ni,ei]
  
def equation_plane(y1, x1, z1, y2, x2, z2, y3, x3, z3):  
  '''Function to find equation of plane.'''
  a1 = x2 - x1 
  b1 = y2 - y1 
  c1 = z2 - z1 
  a2 = x3 - x1 
  b2 = y3 - y1 
  c2 = z3 - z1 
  a = b1 * c2 - b2 * c1 
  b = a2 * c1 - a1 * c2 
  c = a1 * b2 - b1 * a2 
  d = (- a * x1 - b * y1 - c * z1) 
  # plane equation is ax + by + cz + d = 0
  return [a,b,c,d]

def shortest_distance_point_to_plane(y1, x1, z1, a, b, c, d):  
  '''Function to find distance.
     3D point to plane.'''
  d = abs((a * x1 + b * y1 + c * z1 + d))  
  e = math.sqrt(a * a + b * b + c * c)
  return d/e 
    
def distance_point_to_3Dline(y1,x1,z1,y2,x2,z2,y3,x3,z3):
  '''3D point y1,x1,z1
     3D line from y2,x2,z2 to y3,x3,z3'''
  b = math.sqrt((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)
  S = math.sqrt(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))**2 +
                ((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))**2 +
                ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))**2) / 2.0
  return 2*S/b

def calcExtend3DLineCoords(n1,e1,z1,n2,e2,z2,dist):
  lenAB = math.sqrt((n2-n1)**2+(e2-e1)**2+(z2-z1)**2)
  n3=n2+(n2-n1)/lenAB*dist
  e3=e2+(e2-e1)/lenAB*dist
  z3=z2+(z2-z1)/lenAB*dist
  return [n3,e3,z3]

def isect_line_plane_v3_4d(p0, p1, plane, epsilon=1e-6):
  u = sub_v3v3(p1, p0)
  dot = dot_v3v3(plane, u)
  if abs(dot) > epsilon:
    # Calculate a point on the plane
    # (divide can be omitted for unit hessian-normal form).
    p_co=mul_v3_fl(plane,-plane[3]/len_squared_v3(plane))
    w = sub_v3v3(p0, p_co)
    fac = -dot_v3v3(plane, w)/dot
    ang=math.acos(dot/(math.sqrt(len_squared_v3(plane))*math.sqrt(len_squared_v3(u))))
    if ang>math.pi/2.0:
      ang-=math.pi/2.0
    u = mul_v3_fl(u, fac)
    return add_v3v3(p0, u),ang
  else:
    return None,None

def area_by_shoelace(x, y):
  '''Assumes x,y points go around the polygon in one direction'''
  return abs(sum(i * j for i, j in zip(x, y[1:] + y[:1]))
            -sum(i * j for i, j in zip(x[1:] + x[:1], y))) / 2

# generic math functions
def add_v3v3(v0, v1):
  return (
  v0[0] + v1[0],
  v0[1] + v1[1],
  v0[2] + v1[2],
  )

def sub_v3v3(v0, v1):
  return (
    v0[0] - v1[0],
    v0[1] - v1[1],
    v0[2] - v1[2],
  )

def dot_v3v3(v0, v1):
  return (
    (v0[0] * v1[0]) +
    (v0[1] * v1[1]) +
    (v0[2] * v1[2])
    )

def len_squared_v3(v0):
  return dot_v3v3(v0, v0)

def mul_v3_fl(v0, f):
  return (
    v0[0] * f,
    v0[1] * f,
    v0[2] * f,
    )

2 thoughts on “(ฟรี)โปรแกรมภาษาไพทอนบนเครื่องคิดเลขคาสิโอ fx-9750GIII fx-9860GIII และ fx-cg50 โปรแกรมพื้นฐานงานสำรวจชุดที่ 4 (COGO Selected Serie 4)”

Leave a Reply

Your email address will not be published. Required fields are marked *