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

ไม่กี่วันนี้ผมได้ถอยเครื่องคิดเลข fx-9860GIII มาหนึ่งเครื่องราคาประมาณสี่พันห้าร้อยบาท ส่วนน้องๆในที่ทำงานถอย fx-9750GIII มาหนึ่งเครื่องเช่นเดียวกันแต่ราคาย่อมเยากว่า ราคาเครื่องประมาณสามพันบาท สองรุ่นนี้เขียนภาษาไพทอนได้ ไพทอนที่ลงในเครื่องคิดเลขเป็นไพทอนรุ่นเล็กเรียกว่า ไมโครไพทอน (Micropython) แต่ไมโครไพทอนที่ลงในเครื่องคิดเลข ทางคาสิโอลงไลบรารีมาให้ใช้แค่สองไลบรารีคือ math และ random ที่อยากได้มากคือไลบรารี io ที่สามารถเขียนอ่านไฟล์ได้กลับไม่ลงมาให้ ทำให้การใช้งานจำกัดจำเขี่ยเหมือนโดนมัดมือมัดเท้า เอาละไม่เป็นไร ผมจะลองเขียนโปรแกรมไพทอนและพยายามใช้ทรัพยากรอันจำกัดจำเขี่ยนี้ให้ออกมาดูดีที่สุด นี่เป็นความท้าทายเล็กๆน้อยๆอย่างหนึ่ง

ณ ตอนนี้มีเครื่องคิดเลขคาสิโอสามรุ่นที่สนับสนุนภาษาไพทอนเรียงกันดังนี้ fx-9750GIII, fx-9860GIII และ fx-cg50

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

ทดสอบโปรแกรมงานสำรวจที่พัฒนาด้วยภาษาซีบน fx-9750GIII และ fx-9860GIII

ขอย้อนไปโปรแกรมงานสำรวจพื้นฐานบนเครื่องคิดเลข fx-9860GII ที่ผมเขียนด้วยภาษาซี ผมลองเอามาลงในเครื่อง fx-9750III และ fx-9860GIII รันได้ปกติแสดงว่าสถาปัตยกรรมเดียวกัน fx-9860GII ใช้ซีพียู SH4a ส่วนรุ่นอื่นไม่ชัดเจน ความต่างกันที่แน่ๆคือหน่วยความจำ

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

โปรแกรมที่จะเขียนจัดอยู่ในชุด COGO (Coordinate Geometry) ที่มาแก้ปัญหาพิกัดเรขาคณิตในงานสำรวจนั่นเอง ผมเคยเขียนโปรแกรมชุดนี้ด้วยภาษาซีบนเครื่องคิดเลข fx-9860GII (SD) รุ่นก่อน ลองมาดูในภาคภาษาไพทอนบ้าง ผมจะลงซอร์สโค้ดให้ท้ายๆบทความ เผื่อผู้อ่านมีเครื่องคิดเลขรุ่นเหล่านี้จะสามารถก็อปปี้ไปลงเครื่องคิดเลขตัวเองได้

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

สำหรับโปรแกรมพื้นฐานงานสำรวจในชุดนี้จะจัดโปรแกรมย่อยเล็กๆ ไว้ 4 โปรแกรม

  1. Bearing-Dist (2 pt) เมื่อกำหนดจุดค่าพิกัด 2 จุด สำหรับคำนวณหามุมอะซิมัทและระยะทาง
  2. Bearing-Dist(3 pt) เมื่อกำหนดจุดค่าพิกัด 3 จุด สำหรับคำนวณหาง่ามมุมราบ อะซิมัทและระยะทาง ในงานสำรวจก็ได้แก่การตั้งเป้าหลัง  (backsight)  จุดตั้งกล้อง (station) และเป้าหน้า (target)
  3. Coordinate 2D เมื่อกำหนดจุดค่าพิกัด 2 จุด กำหนดมุมราบและระยะราบ คำนวณหาค่าพิกัดจุดที่ 3 คำนวณหาพิกัดจุดที่ 3 การคำนวณคำนวณในระนาบสองมิติอย่างเดียว ไม่มีค่าระดับมาเกี่ยวข้อง
  4. Coordinate 3D เมื่อกำหนดจุดค่าพิกัด 2 จุด กำหนดมุมราบและมุมดิ่ง ระยะทางแบบ slope distance ต้องการคำนวณหาค่าพิกัดและค่าระดับจุดที่ 3

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

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

เลื่อนแถบเคอร์เซอร์ไปที่ COGOSSE1.py กดคีย์ “F1” หรือกดคีย์ “EXE” จะเห็นเมนู

เมนูโปรแกรม พื้นฐานงานสำรวจชุดที่ 1 (COGO Selected Serie 1)

Bearing-Dist (2 pt)

ที่เมนูกดคีย์ “1” และกดคีย์ “EXE” เป็นการคำนวณหาค่ามุมอะซิมัทและระยะทางเมื่อกำหนดจุดค่าพิกัดให้สองจุด ลองทดสอบจากตัวอย่างดังรูป การประยุกต์ใช้งานส่วนใหญ่เป็นตอนที่ช่างสำรวจตั้งกล้องที่หมุดและส่องไปเป้าหลังหรือเป้าหน้าแล้ววัดระยะทางเพื่อตรวจสอบค่าพิกัดเทียบกับค่าที่คำนวณด้วยเครื่องคิดเลข ทดสอบป้อนค่าและจะได้ผลลัพธ์ออกมา จะได้มุมอะซิมัทจากจุดตั้งกล้องไปเป้าหลัง 15°3′ 34.81″ (ตั้งทิศเหนือที่จุดตั้งต้น ที่นี่หมายถึงจุดตั้งกล้องกวาดมุมตามเข็มนาฬิกาไปหาแขนต้นทาง-ปลายทาง) ระยะทาง 87.395 เมตร กดคีย์ “EXE” อีกครั้งโปรแกรมจะออกมาที่เมนูเดิม

แผนผังแสดงตัวอย่างของ Bearing – Dist (2 pt)

Bearing-Dist (3 pt)

ที่เมนูกดคีย์เลข “2” และกดคีย์ “EXE” การประยุกต์ใช้งานส่วนใหญ่จะเป็นการตั้งกล้องส่องไปหมุดเป้าหลังแล้วป้อนค่าพิกัดเป้าหน้าเพื่อตรวจสอบมุมราบหรือไม่ก็จะเป็นการวางผังโดยการเปิดมุมราบและวัดระยะทางที่เป้าบน pole ลองดูตัวอย่างทดสอบ จะได้ผลลัพธ์มาดังนี้ ครั้งแรกจะแสดงมุมอะซิมัทและระยะทางไปเป้าหลังก่อน ถัดไปจะเป็นมุมราบ มุมอะซิมัทและระยะทางไปเป้าหน้า

แผนผังแสดงตัวอย่าง Bearing – Dist (3 pt)

Coordinate 2D

ที่เมนูกดคีย์เลข “3” และกดคีย์ “EXE” เป็นการคำนวณหาค่าพิกัดเป้าหน้าเมื่อกำหนดค่าพิกัดจุดตั้งกล้องและเป้าหลัง กำหนดมุมราบและระยะทาง การคำนวณจะไม่มีค่าระดับมาเกี่ยวข้อง จึงเรียกว่า 2D หรือสองมิติ สำหรับโปรแกรมนี้ผมได้นำค่าสเกลแฟคเตอร์เข้ามาช่วยประยุกต์ใช้ด้วย ในกรณีที่ไม่ต้องใช้ก็ป้อนค่าสเกลแฟคเตอร์นี้ เป็น 1.0

สเกลแฟคเตอร์ตัวนี้แล้วจริงๆคือ Combine Scale Factor (CSF) ที่ได้จาก Elevation Scale Factor (ESF) x Grid Scale Factor (GSF) การประยุกต์ใช้สเกลแฟคเตอร์ส่วนใหญ่นำมาใช้โครงการที่ระบบพิกัดฉากกริดยูทีเอ็มในงานใหญ่ๆยาวๆ เช่นโครงการก่อสร้างถนน รถไฟ เพราะว่าแบบ drawing เราอยู่บนระนาบพิกัดฉาก ให้คิดเสียว่าแบบอยู่บนกระดาษขนาดใหญ่มาตราส่วน 1:1 แล้วเราวัดระยะทางบนผิวโลกที่มีความโค้ง ดังนั้นการวัดระยะทางจะต้องมีการทอนจากบนผิวโค้งเพื่อให้ลงมาเข้ากับระนาบพิกัดฉากของกระดาษ มาลองทดสอบข้อมูล ป้อนข้อมูลค่าพิกัดเป้าหลัง ค่าพิกัดจุดต้องกลอง จากนั้นป้อนมุมป้อนระยะทางและสเกลแฟคเตอร์ (ถ้าใช้) แล้วจะได้ผลลัพธ์เป็นค่าพิกัด N=2642985.886, E=231227.370

แผนผังแสดงตัวอย่าง Coordinate 2D

Coordinate 3D

ที่เมนูกดคีย์ “4” และกดคีย์ “EXE” โปรแกรมคล้าย Coordinate 2D แต่จะมีมิติทางดิ่งเข้ามาเพิ่มดังนั้นที่จุดตั้งกล้องจะวัดความสูงของกล้อง (HI – Height of instrument) และที่เป้าหน้าก็จะต้องวัดความสูงมาด้วย (HT – Height of target) นอกจากนั้นจะมีมุมดิ่ง (Vertical angle) มาด้วย มาดูข้อมูลทดสอบกัน เริ่มจากป้อนค่าพิกัดเป้าหลัง ต่อมาป้อนค่าพิกัดจุดตั้งกล้อง ค่าระดับจุดตั้งกล้อง ความสูงกล้อง ต่อไปป้อนมุมราบ(H.Ang) มุมดิ่ง(V.Ang) ระยะทาง (Slope distance) และความสูงเป้า(HT) และค่าสเกลแฟคเตอร์ (Scale Factor)

โปรแกรมจะคำนวณอะซิมัท ระยะทางจากจุดตั้งกล้องไปเป้าหลัง กดคีย์ “EXE” จะได้ผลลัพธ์ อะซิมัท ระยะราบทั้งระยะบนพื้นโลกและระยะกริดจากจุดตั้งกล้องไปเป้าหน้า สุดท้ายคือผลลัพธ์ที่ต้องการคือค่าพิกัด N=2632215.761, E=232911.782 และค่าระดับ = 6.905 เมตร ของเป้าหน้า

ไดอะแกรมแสดงการวัดหาค่าพิกัดและค่าระดับจากตัวอย่าง

สรุป

ผมเขียนโปรแกรมพื้นฐานงานสำรวจไว้หลายชุด จะทยอยมาลงเป็นซีรี่ย์ ผมเห็นว่ารุ่น fx-9750GIII ราคาประมาณสามพันบาทนั้นไม่ถือว่าแพงมาก สำหรับช่างสำรวจหรือช่างโยธาที่มีการงานเป็นหลักแหล่งสามารถหาซื้อมาใช้งานได้ น้ำหนักเบาจับถนัดมือ

ช่องทางการโปรแกรมมิ่งมีสามช่องทาง ช่องทางที่ 1 คือเขียนด้วยภาษาเบสิคคาสิโอ ก็คล้ายๆกับเครื่องคิดเลขรุ่น fx-5800 อาจจะต่างกันนิดหน่อย ช่องทางที่ 2 คือเขียนด้วยภาษาซี ต้องลงเครื่องมือช่วยพัฒนาโปรแกรม SDK ของคาสิโอ ช่องทางนี้ยากสุดครับ ช่องทางที่ 3 คือเขียนด้วยภาษาไพทอน ที่ผมนำเสนอในบทความนี้ ซึ่งก็ไม่ยากใครมีพื้นฐานโปรแกรมมิ่งภาษาไพทอนนิดๆหน่อยๆก็เขียนได้ ตอนหน้ามาต่อกับโปรแกรมพื้นฐานงานสำรวจชุดที่ 2 (COGO Selected Serie 2) ว่าด้วยเรื่องวงกลมล้วนๆ

มีข้อเตือนใจสำหรับช่างสำรวจสักนิด ปัจจุบันโปรแกรมบนกล้องประมวลผลรวม (Total station) ก็เก่งขึ้นเรื่อยๆ ไม่ต้องกดเครื่องคิดเลขเหมือนสมัยแต่ก่อน ยังมีแอพบนโทรศัพท์มือถือมาช่วยด้วย แต่พื้นฐานจริงๆคือความเข้าใจในเรื่องเรขาคณิต (Geometry) ค่าพิกัด ระยะทาง มุมอะซิมัท จะเกี่ยวข้องกันอย่างแยกกันไม่ออก ดังนั้นต้องทำความเข้าใจหลักวิชาสำรวจให้ลึกซึ้ง เมื่อเข้าใจดีแล้ว การใช้โปรแกรมบนเครื่องคิดเลขหรือบนแอพโทรศัพท์มือถือก็เป็นแค่เครื่องมือที่มาช่วยเหลือให้คำนวณผลลัพธ์มาได้เร็วและถูกต้อง

อนาคตอันใกล้หุ่นยนต์ ปัญญาประดิษฐ์จะมาเร็วกว่าที่คิด หุ่นยนต์หรืออะไรที่มันทำงานแทนเราได้อัตโนมัติ ตัวอย่างคนทำแผนที่ทางอากาศจากโดรนไม่ต้องรู้เรื่อง photogrammetry เลยก็ทำได้ แต่ถ้าคนเราเองไม่ช่วยตัวเองศึกษาค้นคว้า บอกได้คำเดียวว่าเสร็จหุ่นยนต์แน่นอน

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

ผมจะปล่อยซอร์สโค้ดที่อัพเดทไว้ที่เมนูดาวน์โหลด (Download) เมื่อดาวน์โหลดไฟล์ “COGOSSE1.py” และไฟล์ “pbrutils.py” แล้ว สมมุติโหลดมาไว้ที่เครื่องคอมพิวเตอร์ จากนั้นนำเครื่องคิดเลขคาสิโอ fx-9750GIII หรือ fx-9860GIII หรือ fx-cg50 มาเสียบเข้าสาย USB ต่อกับเครื่องคอมพิวเตอร์ ให้กดคีย์ “F1” เพื่อจำลองเป็นไดรว์ จากนั้น copy สองไฟล์นี้เข้าไปที่ไดรว์ของเครื่องคิดเลข รายละเอียดอ่านที่ลิ๊งค์ที่ผมจัดทำมาเป็นพิเศษ

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

from pbrutils import split_angle,deg2dms,DMS2str
from pbrutils import calcAzimuth,calcDistance,calcDiff2Angle
from pbrutils import calcNextAzimuth,calcNextCoordinate2D
from pbrutils import calcNextCoordinate3D
from pbrutils import PI,DEG2RAD,RAD2DEG,clear
import math
 
def print_menu():
  print('COGO Selected Serie 1')
  print('1: Bearing-Dist(2pt)')
  print('2: Bearing-Dist(3pt)')
  print('3: Coordinate 2D')
  print('4: Coordinate 3D')
  print('5: Exit')

def bearing_dist_2pt():
  nbs=float(input("N-BS=" ))
  ebs=float(input("E-BS=" ))
  nin=float(input("N-Instr=" ))
  ein=float(input("E-Instr=" ))
  azi21=calcAzimuth(nin,ein,nbs,ebs)
  dist21=calcDistance(nin,ein,nbs,ebs)
  dazi21=azi21*RAD2DEG
  print(6*'-',"Output",6*'-')
  print("Instrument => BS")
  sn,d,m,s=deg2dms(dazi21)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st) 
  print("Dist=","{0:.3f} m".format(dist21)) 
  input("Press EXE key")  

def bearing_dist_3pt():
  nbs=float(input("N-BS=" ))
  ebs=float(input("E-BS=" ))
  nin=float(input("N-Instr=" ))
  ein=float(input("E-Instr=" ))
  nfs=float(input("N-FS=" ))
  efs=float(input("E-FS=" ))
  azi21=calcAzimuth(nin,ein,nbs,ebs)
  dist21=calcDistance(nin,ein,nbs,ebs)
  azi23=calcAzimuth(nin,ein,nfs,efs)
  dist23=calcDistance(nin,ein,nfs,efs)
  dazi21=azi21*RAD2DEG
  dazi23=azi23*RAD2DEG
  print(6*'-',"Output",6*'-')
  print("Instrument => BS")
  sn,d,m,s=deg2dms(dazi21)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st) 
  print("Dist=","{0:.3f} m".format(dist21)) 
  input("Press EXE key")
  print("Instrument => FS")
  sn,d,m,s=deg2dms(dazi23)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st) 
  print("Dist=","{0:.3f} m".format(dist23)) 
  hang=RAD2DEG*calcDiff2Angle(azi23,azi21)
  sn,d,m,s=deg2dms(hang)
  st=DMS2str(d,m,s,2)
  print("H Angle=",st)  
  input("Press EXE key")  
  
def cooordinate2D():
  nbs=float(input("N-BS=?"))
  ebs=float(input("E-BS=?"))
  nin=float(input("N-Instr=?"))
  ein=float(input("E-Instr=?"))
  ang=input("H Angle=?")
  hang=split_angle(ang)*DEG2RAD
  grndist=float(input("Ground Dist=?"))
  sf=float(input("SF=?"))
  grddist=sf*grndist
  azi12=calcAzimuth(nbs,ebs,nin,ein)
  azi21=calcAzimuth(nin,ein,nbs,ebs)
  dist21=calcDistance(nin,ein,nbs,ebs)  
  nt,et=calcNextCoordinate2D(nin,ein,azi12,hang,grddist)
  azi23=calcNextAzimuth(hang,azi12) 
  print(6*'-',"Output",6*'-') 
  print("Instrument => BS")
  sn,d,m,s=deg2dms(azi21*RAD2DEG)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st)
  print("Grd Dist=","{0:.3f} m".format(dist21)) 
  print("Instrument => FS")
  sn,d,m,s=deg2dms(azi23*RAD2DEG)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st) 
  print("Gnd Dist=","{0:.3f} m".format(grndist)) 
  print("Grd Dist=","{0:.3f} m".format(grddist)) 
  input("Press EXE key") 
  print("N=","{0:.3f} m".format(nt))
  print("E=","{0:.3f} m".format(et))
  input("Press EXE key") 

def cooordinate3D():
  nbs=float(input("N-BS=?"))
  ebs=float(input("E-BS=?"))
  nin=float(input("N-Instr=?"))
  ein=float(input("E-Instr=?"))
  ele=float(input("Elev=?"))
  hi=float(input("HI=?"))
  ang=input("H Angle=?")
  hang=split_angle(ang)*DEG2RAD
  ang=input("V Angle=?")
  vang=split_angle(ang)*DEG2RAD
  sgrndist=float(input("Slope Dist=?"))
  sf=1.0
  sf=float(input("SF=?"))
  if sf<1.0:
    sf=1.0
  sgrddist=sgrndist*sf
  ht=float(input("HT=?"))
  azi21=calcAzimuth(nin,ein,nbs,ebs)
  azi12=calcAzimuth(nbs,ebs,nin,ein)
  dist21=calcDistance(nin,ein,nbs,ebs)
  azi23=calcNextAzimuth(hang,azi12)
  nt,et,zt=calcNextCoordinate3D(nin,ein,ele,azi12,hang,vang,sgrddist,hi,ht)
  print(6*'-',"Output",6*'-') 
  print("Instrument => BS")
  sn,d,m,s=deg2dms(azi21*RAD2DEG)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st)
  print("Grd Dist=","{0:.3f} m".format(dist21)) 
  print("Instrument => FS")
  sn,d,m,s=deg2dms(azi23*RAD2DEG)
  st=DMS2str(d,m,s,2)
  print("Azimuth=",st) 
  input("Press EXE key") 
  print("S Gnd Dist= {0:.3f} m".format(sgrndist)) 
  print("S Grd Dist= {0:.3f} m".format(sgrddist)) 
  print("H Gnd Dist= {:.3f} m".format(sgrndist*math.sin(vang)))
  print("H Grd Dist= {:.3f} m".format(sgrddist*math.sin(vang)))
  print("Height diff= {:.3f} m".format(zt-ele))
  input("Press EXE key")
  print("N=","{0:.3f} m".format(nt))
  print("E=","{0:.3f} m".format(et))  
  print("Elev=","{0:.3f} m".format(zt)) 
  input("Press EXE key") 
  
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
    bearing_dist_2pt()
    ex_choice=1
  elif (choice==2):
    loop=True
    bearing_dist_3pt()
    ex_choice=2
  elif (choice==3):
    cooordinate2D()
    loop=True
    ex_choice=3
  elif (choice==4):
    loop=True
    cooordinate3D()
    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,
    )

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

Leave a Reply

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