github.com/archlabjp/eeslism-go@v0.0.0-20231109122333-4bb7bfcdf292/eeslism/COORDNT.go

//This file is part of EESLISM.
//
//Foobar is free software : you can redistribute itand /or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//Foobar is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the
//GNU General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with Foobar.If not, see < https://www.gnu.org/licenses/>.

/*

	入力データの計算用構造体への変換
	FILE=COORDNT.c
	Create Date=1999.6.7
*/

package eeslism

import (
	"fmt"
	"math"
	"os"
)

func LP_COORDNT(
	lpn *int,
	bdpn int,
	obsn int,
	treen int,
	polyn int,
	poly []POLYGN,
	tree []TREE,
	obs []OBS,
	BDP []BBDP,
	lp []P_MENN,
) {
	var name string
	var i, j, k, p int
	var sWa, cWa, sWb, cWb, a, b, c, chWA, shWA float64
	var csWA, ssWA, csWb, ssWb float64
	var caWb, saWa, caWa, saWb float64
	var cbWa, sbWa float64

	k = 0
	const M_rad = math.Pi / 180.0
	for i = 0; i < bdpn; i++ {
		sWb = math.Sin(BDP[i].Wb * M_rad)
		cWb = math.Cos(BDP[i].Wb * M_rad)
		sWa = math.Sin(-BDP[i].Wa * M_rad)
		cWa = math.Cos(-BDP[i].Wa * M_rad)

		if BDP[i].sumsblk != 0 {
			for j = 0; j < BDP[i].sumsblk; j++ {

				if BDP[i].SBLK[j].sbfname == "HISASI" {
					a = 0.0
					b = 0.0
					c = 0.0
					chWA = 0.0
					shWA = 0.0

					a = BDP[i].x0 + BDP[i].SBLK[j].x*cWa - BDP[i].SBLK[j].y*cWb*sWa
					b = BDP[i].y0 + BDP[i].SBLK[j].x*sWa + BDP[i].SBLK[j].y*cWb*cWa
					c = BDP[i].z0 + BDP[i].SBLK[j].y*sWb

					chWA = math.Cos((BDP[i].SBLK[j].WA - BDP[i].Wb) * M_rad)
					shWA = math.Sin((BDP[i].SBLK[j].WA - BDP[i].Wb) * M_rad)

					lp[k].opname = BDP[i].SBLK[j].snbname
					lp[k].wa = BDP[i].Wa
					lp[k].wb = BDP[i].Wb + (180.0 - BDP[i].SBLK[j].WA)
					if lp[k].wb > 180.0 {
						lp[k].wb = 360.0 - lp[k].wb
					}

					lp[k].wd = 0
					lp[k].ref = 0.0 /*--付設障害物からの反射率を０としている-*/
					for p = 0; p < 366; p++ {
						lp[k].shad[p] = 1
					}

					lp[k].sbflg = 1

					lp[k].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k].polyd = 4
					lp[k].P = make([]XYZ, 4)

					lp[k].P[0].X = a
					lp[k].P[0].Y = b
					lp[k].P[0].Z = c

					lp[k].P[1].X = a + BDP[i].SBLK[j].D*chWA*sWa
					lp[k].P[1].Y = b - BDP[i].SBLK[j].D*chWA*cWa
					lp[k].P[1].Z = c + BDP[i].SBLK[j].D*shWA

					lp[k].P[2].X = lp[k].P[1].X + BDP[i].SBLK[j].W*cWa
					lp[k].P[2].Y = lp[k].P[1].Y + BDP[i].SBLK[j].W*sWa
					lp[k].P[2].Z = lp[k].P[1].Z

					lp[k].P[3].X = a + BDP[i].SBLK[j].W*cWa
					lp[k].P[3].Y = b + BDP[i].SBLK[j].W*sWa
					lp[k].P[3].Z = c

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
					lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
					lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
					CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

					//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

					k = k + 1
				} else if BDP[i].SBLK[j].sbfname == "BARUKONI" {

					a = 0.0
					b = 0.0
					c = 0.0
					chWA = 0.0
					shWA = 0.0

					a = BDP[i].x0 + BDP[i].SBLK[j].x*cWa - BDP[i].SBLK[j].y*cWb*sWa
					b = BDP[i].y0 + BDP[i].SBLK[j].x*sWa + BDP[i].SBLK[j].y*cWb*cWa
					c = BDP[i].z0 + BDP[i].SBLK[j].y*sWb

					chWA = math.Cos((90.0 - BDP[i].Wb) * M_rad)
					shWA = math.Sin((90.0 - BDP[i].Wb) * M_rad)

					lp[k].opname = BDP[i].SBLK[j].snbname

					lp[k].wa = BDP[i].Wa
					lp[k].wb = BDP[i].Wb + 90.0
					if lp[k].wb > 180.0 {
						lp[k].wb = -(360.0 - lp[k].wb)
					}

					lp[k].ref = 0.0
					lp[k].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k].shad[p] = 1
					}

					lp[k].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k].polyd = 4
					lp[k].P = make([]XYZ, 4)

					lp[k].wd = 0

					lp[k].P[0].X = a
					lp[k].P[0].Y = b
					lp[k].P[0].Z = c

					lp[k].P[1].X = a + BDP[i].SBLK[j].D*chWA*sWa
					lp[k].P[1].Y = b - BDP[i].SBLK[j].D*chWA*cWa
					lp[k].P[1].Z = c + BDP[i].SBLK[j].D*shWA

					lp[k].P[2].X = lp[k].P[1].X + BDP[i].SBLK[j].W*cWa
					lp[k].P[2].Y = lp[k].P[1].Y + BDP[i].SBLK[j].W*sWa
					lp[k].P[2].Z = lp[k].P[1].Z

					lp[k].P[3].X = a + BDP[i].SBLK[j].W*cWa
					lp[k].P[3].Y = b + BDP[i].SBLK[j].W*sWa
					lp[k].P[3].Z = c

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
					lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
					lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
					CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

					//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

					name = "2" + lp[k].opname
					lp[k+1].opname = name

					lp[k+1].wa = BDP[i].Wa - 90.0
					if lp[k+1].wa <= -180.0 {
						lp[k+1].wa = 360.0 + lp[k+1].wa
					}
					lp[k+1].wb = 90.0
					lp[k+1].ref = 0.0
					lp[k+1].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k+1].shad[p] = 1
					}

					lp[k+1].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k+1].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k+1].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k+1].polyd = 4
					lp[k+1].P = make([]XYZ, 4)

					lp[k+1].wd = 0

					lp[k+1].P[0] = lp[k].P[0]
					lp[k+1].P[1] = lp[k].P[1]
					lp[k+1].P[2].X = lp[k+1].P[1].X + BDP[i].SBLK[j].H*cWb*sWa
					lp[k+1].P[2].Y = lp[k+1].P[1].Y - BDP[i].SBLK[j].H*cWb*cWa
					lp[k+1].P[2].Z = lp[k+1].P[1].Z - BDP[i].SBLK[j].H*sWb
					lp[k+1].P[3].X = a + BDP[i].SBLK[j].H*cWb*sWa
					lp[k+1].P[3].Y = b - BDP[i].SBLK[j].H*cWb*cWa
					lp[k+1].P[3].Z = c - BDP[i].SBLK[j].H*sWb

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k+1].e.Z = math.Cos(lp[k+1].wb * M_rad)
					lp[k+1].e.Y = -math.Sin(lp[k+1].wb*M_rad) * math.Cos(lp[k+1].wa*M_rad)
					lp[k+1].e.X = -math.Sin(lp[k+1].wb*M_rad) * math.Sin(lp[k+1].wa*M_rad)
					CAT(&lp[k+1].e.X, &lp[k+1].e.Y, &lp[k+1].e.Z)

					//HOUSEN2(&lp[k+1].P[0],&lp[k+1].P[1],&lp[k+1].P[2],&lp[k+1].e)

					name = "3" + lp[k].opname
					lp[k+2].opname = name

					lp[k+2].wa = BDP[i].Wa
					lp[k+2].wb = BDP[i].Wb - 90.0
					lp[k+2].ref = BDP[i].SBLK[j].ref
					lp[k+2].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k+2].shad[p] = 1
					}

					lp[k+2].wd = 0

					lp[k+2].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k+2].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k+2].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k+2].polyd = 4
					lp[k+2].P = make([]XYZ, 4)

					lp[k+2].P[0] = lp[k+1].P[3]
					lp[k+2].P[1] = lp[k+1].P[2]
					lp[k+2].P[2].X = lp[k+2].P[1].X + BDP[i].SBLK[j].W*cWa
					lp[k+2].P[2].Y = lp[k+2].P[1].Y + BDP[i].SBLK[j].W*sWa
					lp[k+2].P[2].Z = lp[k+2].P[1].Z
					lp[k+2].P[3].X = lp[k+2].P[0].X + BDP[i].SBLK[j].W*cWa
					lp[k+2].P[3].Y = lp[k+2].P[0].Y + BDP[i].SBLK[j].W*sWa
					lp[k+2].P[3].Z = lp[k+2].P[0].Z

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k+2].e.Z = math.Cos(lp[k+2].wb * M_rad)
					lp[k+2].e.Y = -math.Sin(lp[k+2].wb*M_rad) * math.Cos(lp[k+2].wa*M_rad)
					lp[k+2].e.X = -math.Sin(lp[k+2].wb*M_rad) * math.Sin(lp[k+2].wa*M_rad)
					CAT(&lp[k+2].e.X, &lp[k+2].e.Y, &lp[k+2].e.Z)

					//HOUSEN2(&lp[k+2].P[0],&lp[k+2].P[1],&lp[k+2].P[2],&lp[k+2].e) ;

					name = "4" + lp[k].opname
					lp[k+3].opname = name

					lp[k+3].wa = BDP[i].Wa + 90.0
					if lp[k+3].wa >= 180.0 {
						lp[k+3].wa = lp[k+3].wa - 360.0
					}
					lp[k+3].wb = 90.0
					lp[k+3].ref = 0.0
					lp[k+3].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k+3].shad[p] = 1
					}

					lp[k+3].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k+3].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k+3].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k+3].polyd = 4
					lp[k+3].P = make([]XYZ, 4)

					lp[k+3].wd = 0

					lp[k+3].P[0] = lp[k].P[3]
					lp[k+3].P[1] = lp[k].P[2]
					lp[k+3].P[2] = lp[k+2].P[2]
					lp[k+3].P[3] = lp[k+2].P[3]

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k+3].e.Z = math.Cos(lp[k+3].wb * M_rad)
					lp[k+3].e.Y = -math.Sin(lp[k+3].wb*M_rad) * math.Cos(lp[k+3].wa*M_rad)
					lp[k+3].e.X = -math.Sin(lp[k+3].wb*M_rad) * math.Sin(lp[k+3].wa*M_rad)
					CAT(&lp[k+3].e.X, &lp[k+3].e.Y, &lp[k+3].e.Z)

					//HOUSEN2(&lp[k+3].P[0],&lp[k+3].P[1],&lp[k+3].P[2],&lp[k+3].e)

					name = "5" + lp[k].opname
					lp[k+4].opname = name

					lp[k+4].wa = BDP[i].Wa - 180
					if lp[k+4].wa > 180 {
						lp[k+4].wa = lp[k+4].wa - 360
					} else if lp[k+4].wa < 180 {
						lp[k+4].wa = lp[k+4].wa + 360
					}

					lp[k+4].wb = BDP[i].Wb
					lp[k+4].ref = 0.0
					lp[k+4].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k+4].shad[p] = 1
					}

					lp[k+4].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k+4].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k+4].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k+4].polyd = 4
					lp[k+4].P = make([]XYZ, 4)

					lp[k+4].wd = 0

					lp[k+4].P[0] = lp[k+1].P[2]
					lp[k+4].P[1].X = lp[k+4].P[0].X - BDP[i].SBLK[j].h*sWa*cWb
					lp[k+4].P[1].Y = lp[k+4].P[0].Y + BDP[i].SBLK[j].h*cWa*cWb
					lp[k+4].P[1].Z = lp[k+4].P[0].Z + BDP[i].SBLK[j].h*sWb
					lp[k+4].P[2].X = lp[k+4].P[1].X + BDP[i].SBLK[j].W*cWa
					lp[k+4].P[2].Y = lp[k+4].P[1].Y + BDP[i].SBLK[j].W*sWa
					lp[k+4].P[2].Z = lp[k+4].P[1].Z
					lp[k+4].P[3] = lp[k+3].P[2]

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k+4].e.Z = math.Cos(lp[k+4].wb * M_rad)
					lp[k+4].e.Y = -math.Sin(lp[k+4].wb*M_rad) * math.Cos(lp[k+4].wa*M_rad)
					lp[k+4].e.X = -math.Sin(lp[k+4].wb*M_rad) * math.Sin(lp[k+4].wa*M_rad)
					CAT(&lp[k+4].e.X, &lp[k+4].e.Y, &lp[k+4].e.Z)

					//HOUSEN2(&lp[k+4].P[0],&lp[k+4].P[1],&lp[k+4].P[2],&lp[k+4].e)
					k = k + 5
				} else if BDP[i].SBLK[j].sbfname == "SODEKABE" {

					a = 0.0
					b = 0.0
					c = 0.0
					csWA = 0.0
					ssWA = 0.0
					csWb = 0.0
					ssWb = 0.0

					a = BDP[i].x0 + BDP[i].SBLK[j].x*cWa - BDP[i].SBLK[j].y*cWb*sWa
					b = BDP[i].y0 + BDP[i].SBLK[j].x*sWa + BDP[i].SBLK[j].y*cWb*cWa
					c = BDP[i].z0 + BDP[i].SBLK[j].y*sWb
					csWA = math.Cos((-BDP[i].Wa - BDP[i].SBLK[j].WA) * M_rad)
					ssWA = math.Sin((-BDP[i].Wa - BDP[i].SBLK[j].WA) * M_rad)
					csWb = math.Cos((90.0 - BDP[i].Wb) * M_rad)
					ssWb = math.Sin((90.0 - BDP[i].Wb) * M_rad)

					lp[k].opname = BDP[i].SBLK[j].snbname

					lp[k].wa = BDP[i].Wa - 90.0
					if lp[k].wa <= -180.0 {
						lp[k].wa = 360.0 + lp[k].wa
					}
					lp[k].wb = 90.0
					lp[k].ref = 0.0
					lp[k].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k].shad[p] = 1
					}

					lp[k].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k].polyd = 4
					lp[k].P = make([]XYZ, 4)

					lp[k].wd = 0

					lp[k].P[0].X = a
					lp[k].P[0].Y = b
					lp[k].P[0].Z = c
					lp[k].P[1].X = a + BDP[i].SBLK[j].D*csWb*csWA
					lp[k].P[1].Y = b + BDP[i].SBLK[j].D*csWb*ssWA
					lp[k].P[1].Z = c + BDP[i].SBLK[j].D*ssWb
					lp[k].P[2].X = lp[k].P[1].X + BDP[i].SBLK[j].H*cWb*sWa
					lp[k].P[2].Y = lp[k].P[1].Y - BDP[i].SBLK[j].H*cWb*cWa
					lp[k].P[2].Z = lp[k].P[1].Z - BDP[i].SBLK[j].H*sWb
					lp[k].P[3].X = a + BDP[i].SBLK[j].H*cWb*sWa
					lp[k].P[3].Y = b - BDP[i].SBLK[j].H*cWb*cWa
					lp[k].P[3].Z = c - BDP[i].SBLK[j].H*sWb

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
					lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
					lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
					CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

					//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e) ;

					k = k + 1

				} else if BDP[i].SBLK[j].sbfname == "MADOHIYOKE" {

					a = 0.0
					b = 0.0
					c = 0.0
					chWA = 0.0
					shWA = 0.0

					a = BDP[i].x0 + BDP[i].SBLK[j].x*cWa - BDP[i].SBLK[j].y*cWb*sWa
					b = BDP[i].y0 + BDP[i].SBLK[j].x*sWa + BDP[i].SBLK[j].y*cWb*cWa
					c = BDP[i].z0 + BDP[i].SBLK[j].y*sWb
					chWA = math.Cos((90.0 - BDP[i].Wb) * M_rad)
					shWA = math.Sin((90.0 - BDP[i].Wb) * M_rad)

					lp[k].opname = BDP[i].SBLK[j].snbname

					lp[k].wa = BDP[i].Wa - 180
					if lp[k].wa > 180 {
						lp[k].wa = lp[k].wa - 360
					} else if lp[k].wa < 180 {
						lp[k].wa = lp[k].wa + 360
					}

					lp[k].wb = BDP[i].Wb
					lp[k].ref = 0.0
					lp[k].sbflg = 1

					for p = 0; p < 366; p++ {
						lp[k].shad[p] = 1
					}

					lp[k].rgb[0] = BDP[i].SBLK[j].rgb[0]
					lp[k].rgb[1] = BDP[i].SBLK[j].rgb[1]
					lp[k].rgb[2] = BDP[i].SBLK[j].rgb[2]
					lp[k].polyd = 4
					lp[k].P = make([]XYZ, 4)

					lp[k].wd = 0

					lp[k].P[0].X = a + BDP[i].SBLK[j].D*chWA*sWa
					lp[k].P[0].Y = b - BDP[i].SBLK[j].D*chWA*cWa
					lp[k].P[0].Z = c + BDP[i].SBLK[j].H*shWA
					lp[k].P[1].X = lp[k].P[0].X + BDP[i].SBLK[j].H*cWb*sWa
					lp[k].P[1].Y = lp[k].P[0].Y - BDP[i].SBLK[j].H*cWb*cWa
					lp[k].P[1].Z = lp[k].P[0].Z - BDP[i].SBLK[j].H*sWb
					lp[k].P[3].X = lp[k].P[0].X + BDP[i].SBLK[j].W*cWa
					lp[k].P[3].Y = lp[k].P[0].Y + BDP[i].SBLK[j].W*sWa
					lp[k].P[3].Z = lp[k].P[0].Z
					lp[k].P[2].X = lp[k].P[1].X + BDP[i].SBLK[j].W*cWa
					lp[k].P[2].Y = lp[k].P[1].Y + BDP[i].SBLK[j].W*sWa
					lp[k].P[2].Z = lp[k].P[1].Z

					/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
					lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
					lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
					lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
					CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

					//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e) ;
					k = k + 1
				}
			}
		}

	}

	/*-----------------------------------------------------*/

	for i = 0; i < obsn; i++ {
		if obs[i].fname == "rect" {
			caWb = 0.0
			saWa = 0.0
			caWa = 0.0
			saWb = 0.0

			caWb = math.Cos(obs[i].Wb * M_rad)
			saWa = math.Sin(-obs[i].Wa * M_rad)
			saWb = math.Sin(obs[i].Wb * M_rad)
			caWa = math.Cos(-obs[i].Wa * M_rad)

			lp[k].opname = obs[i].obsname

			lp[k].wa = obs[i].Wa
			lp[k].wb = obs[i].Wb
			lp[k].ref = obs[i].ref[0]
			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = obs[i].rgb[0]
			lp[k].rgb[1] = obs[i].rgb[1]
			lp[k].rgb[2] = obs[i].rgb[2]
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].wd = 0
			lp[k].sbflg = 0

			lp[k].P[0].X = obs[i].x
			lp[k].P[0].Y = obs[i].y
			lp[k].P[0].Z = obs[i].z
			lp[k].P[1].X = obs[i].x - obs[i].H*caWb*saWa
			lp[k].P[1].Y = obs[i].y + obs[i].H*caWb*caWa
			lp[k].P[1].Z = obs[i].z + obs[i].H*saWb
			lp[k].P[2].X = obs[i].x + obs[i].W*caWa - obs[i].H*caWb*saWa
			lp[k].P[2].Y = obs[i].y + obs[i].H*caWb*caWa + obs[i].W*saWa
			lp[k].P[2].Z = obs[i].z + obs[i].H*saWb
			lp[k].P[3].X = obs[i].x + obs[i].W*caWa
			lp[k].P[3].Y = obs[i].y + obs[i].W*saWa
			lp[k].P[3].Z = obs[i].z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1
		} else if obs[i].fname == "cube" {
			cbWa = 0.0
			sbWa = 0.0

			cbWa = math.Cos(-obs[i].Wa * M_rad)
			sbWa = math.Sin(-obs[i].Wa * M_rad)

			lp[k].opname = obs[i].obsname

			lp[k].wa = obs[i].Wa
			lp[k].wb = 90.0
			lp[k].ref = obs[i].ref[0]
			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = obs[i].rgb[0]
			lp[k].rgb[1] = obs[i].rgb[1]
			lp[k].rgb[2] = obs[i].rgb[2]
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].wd = 0
			lp[k].sbflg = 0

			lp[k].P[0].X = obs[i].x
			lp[k].P[0].Y = obs[i].y
			lp[k].P[0].Z = obs[i].z
			lp[k].P[1].X = obs[i].x
			lp[k].P[1].Y = obs[i].y
			lp[k].P[1].Z = obs[i].z + obs[i].H
			lp[k].P[2].X = obs[i].x + obs[i].W*cbWa
			lp[k].P[2].Y = obs[i].y + obs[i].W*sbWa
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = lp[k].P[2].X
			lp[k].P[3].Y = lp[k].P[2].Y
			lp[k].P[3].Z = lp[k].P[0].Z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			name = "2" + lp[k].opname
			lp[k+1].opname = name

			lp[k+1].wa = obs[i].Wa - 90.0
			if lp[k+1].wa <= -180.0 {
				lp[k+1].wa = 360.0 + lp[k+1].wa
			}
			lp[k+1].wb = 90.0
			lp[k+1].ref = obs[i].ref[1]
			for p = 0; p < 366; p++ {
				lp[k+1].shad[p] = 1
			}

			lp[k+1].rgb[0] = obs[i].rgb[0]
			lp[k+1].rgb[1] = obs[i].rgb[1]
			lp[k+1].rgb[2] = obs[i].rgb[2]
			lp[k+1].polyd = 4
			lp[k+1].P = make([]XYZ, 4)

			lp[k+1].wd = 0
			lp[k+1].sbflg = 0

			lp[k+1].P[0] = lp[k].P[3]
			lp[k+1].P[1] = lp[k].P[2]
			lp[k+1].P[2].X = lp[k+1].P[1].X - obs[i].D*sbWa
			lp[k+1].P[2].Y = lp[k+1].P[1].Y + obs[i].D*cbWa
			lp[k+1].P[2].Z = lp[k].P[2].Z
			lp[k+1].P[3].X = lp[k+1].P[2].X
			lp[k+1].P[3].Y = lp[k+1].P[2].Y
			lp[k+1].P[3].Z = lp[k+1].P[0].Z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k+1].e.Z = math.Cos(lp[k+1].wb * M_rad)
			lp[k+1].e.Y = -math.Sin(lp[k+1].wb*M_rad) * math.Cos(lp[k+1].wa*M_rad)
			lp[k+1].e.X = -math.Sin(lp[k+1].wb*M_rad) * math.Sin(lp[k+1].wa*M_rad)
			CAT(&lp[k+1].e.X, &lp[k+1].e.Y, &lp[k+1].e.Z)
			//HOUSEN2(&lp[k+1].P[0],&lp[k+1].P[1],&lp[k+1].P[2],&lp[k+1].e)

			name = "3" + lp[k].opname
			lp[k+2].opname = name

			lp[k+2].wa = obs[i].Wa - 180.0
			if lp[k+2].wa <= -180.0 {
				lp[k+2].wa = 360.0 + lp[k+2].wa
			}
			lp[k+2].wb = 90.0
			lp[k+2].ref = obs[i].ref[2]
			for p = 0; p < 366; p++ {
				lp[k+2].shad[p] = 1
			}

			lp[k+2].rgb[0] = obs[i].rgb[0]
			lp[k+2].rgb[1] = obs[i].rgb[1]
			lp[k+2].rgb[2] = obs[i].rgb[2]
			lp[k+2].polyd = 4
			lp[k+2].P = make([]XYZ, 4)

			lp[k+2].wd = 0
			lp[k+2].sbflg = 0
			lp[k+2].P[0] = lp[k+1].P[3]
			lp[k+2].P[1] = lp[k+1].P[2]
			lp[k+2].P[2].X = obs[i].x - obs[i].D*sbWa
			lp[k+2].P[2].Y = obs[i].y + obs[i].D*cbWa
			lp[k+2].P[2].Z = lp[k+2].P[1].Z
			lp[k+2].P[3].X = lp[k+2].P[2].X
			lp[k+2].P[3].Y = lp[k+2].P[2].Y
			lp[k+2].P[3].Z = lp[k+2].P[0].Z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k+2].e.Z = math.Cos(lp[k+2].wb * M_rad)
			lp[k+2].e.Y = -math.Sin(lp[k+2].wb*M_rad) * math.Cos(lp[k+2].wa*M_rad)
			lp[k+2].e.X = -math.Sin(lp[k+2].wb*M_rad) * math.Sin(lp[k+2].wa*M_rad)
			CAT(&lp[k+2].e.X, &lp[k+2].e.Y, &lp[k+2].e.Z)
			//HOUSEN2(&lp[k+2].P[0],&lp[k+2].P[1],&lp[k+2].P[2],&lp[k+2].e)

			name = "4" + lp[k].opname
			lp[k+3].opname = name

			lp[k+3].wa = obs[i].Wa + 90.0
			if lp[k+3].wa > 180.0 {
				lp[k+3].wa = 360.0 - lp[k+3].wa
			}
			lp[k+3].wb = 90.0
			lp[k+3].ref = obs[i].ref[3]
			for p = 0; p < 366; p++ {
				lp[k+3].shad[p] = 1
			}

			lp[k+3].rgb[0] = obs[i].rgb[0]
			lp[k+3].rgb[1] = obs[i].rgb[1]
			lp[k+3].rgb[2] = obs[i].rgb[2]
			lp[k+3].polyd = 4
			lp[k+3].P = make([]XYZ, 4)

			lp[k+3].wd = 0
			lp[k+3].sbflg = 0
			lp[k+3].P[0] = lp[k+2].P[3]
			lp[k+3].P[1] = lp[k+2].P[2]
			lp[k+3].P[2] = lp[k].P[1]
			lp[k+3].P[3] = lp[k].P[0]

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k+3].e.Z = math.Cos(lp[k+3].wb * M_rad)
			lp[k+3].e.Y = -math.Sin(lp[k+3].wb*M_rad) * math.Cos(lp[k+3].wa*M_rad)
			lp[k+3].e.X = -math.Sin(lp[k+3].wb*M_rad) * math.Sin(lp[k+3].wa*M_rad)
			CAT(&lp[k+3].e.X, &lp[k+3].e.Y, &lp[k+3].e.Z)
			//HOUSEN2(&lp[k+3].P[0],&lp[k+3].P[1],&lp[k+3].P[2],&lp[k+3].e)

			k = k + 4
		} else if obs[i].fname == "r_tri" {

			caWb = 0.0
			saWa = 0.0
			caWa = 0.0
			saWb = 0.0

			caWb = math.Cos(obs[i].Wb * M_rad)
			saWa = math.Sin(obs[i].Wa * M_rad)
			saWb = math.Sin(obs[i].Wb * M_rad)
			caWa = math.Cos(obs[i].Wa * M_rad)

			lp[k].opname = obs[i].obsname

			lp[k].wa = obs[i].Wa
			lp[k].wb = obs[i].Wb
			lp[k].ref = obs[i].ref[0]
			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = obs[i].rgb[0]
			lp[k].rgb[1] = obs[i].rgb[1]
			lp[k].rgb[2] = obs[i].rgb[2]
			lp[k].polyd = 3
			lp[k].P = make([]XYZ, 3)

			lp[k].wd = 0
			lp[k].sbflg = 0

			lp[k].P[0].X = obs[i].x
			lp[k].P[0].Y = obs[i].y
			lp[k].P[0].Z = obs[i].z
			lp[k].P[1].X = obs[i].x - obs[i].H*caWb*saWa
			lp[k].P[1].Y = obs[i].y + obs[i].H*caWb*caWa
			lp[k].P[1].Z = obs[i].z + obs[i].H*saWb
			lp[k].P[2].X = obs[i].x + obs[i].W*caWa
			lp[k].P[2].Y = obs[i].y + obs[i].W*saWa
			lp[k].P[2].Z = obs[i].z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1
		} else if obs[i].fname == "i_tri" {

			caWb = 0.0
			saWa = 0.0
			caWa = 0.0
			saWb = 0.0

			caWb = math.Cos(obs[i].Wb * M_rad)
			saWa = math.Sin(-obs[i].Wa * M_rad)
			saWb = math.Sin(obs[i].Wb * M_rad)
			caWa = math.Cos(-obs[i].Wa * M_rad)

			lp[k].opname = obs[i].obsname

			lp[k].wa = obs[i].Wa
			lp[k].wb = obs[i].Wb
			lp[k].ref = obs[i].ref[0]
			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = obs[i].rgb[0]
			lp[k].rgb[1] = obs[i].rgb[1]
			lp[k].rgb[2] = obs[i].rgb[2]
			lp[k].polyd = 3
			lp[k].P = make([]XYZ, 3)

			lp[k].wd = 0
			lp[k].sbflg = 0

			lp[k].P[0].X = obs[i].x
			lp[k].P[0].Y = obs[i].y
			lp[k].P[0].Z = obs[i].z
			lp[k].P[1].X = obs[i].x + ((obs[i].W)/2)*caWa - obs[i].H*caWb*saWa
			lp[k].P[1].Y = obs[i].y + obs[i].H*caWb*caWa + ((obs[i].W)/2)*saWa
			lp[k].P[1].Z = obs[i].z + obs[i].H*saWb
			lp[k].P[2].X = obs[i].x + obs[i].W*caWa
			lp[k].P[2].Y = obs[i].y + obs[i].W*saWa
			lp[k].P[2].Z = obs[i].z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1
		} else {
			fmt.Printf("error--**COORDNT-lp\n")
			os.Exit(1)
		}

	}

	/*--------------------------------------------------------*/
	for i = 0; i < treen; i++ {
		if tree[i].treetype == "treeA" {
			name = tree[i].treename + "-m1"
			lp[k].opname = name

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.4
			lp[k].rgb[1] = 0.3
			lp[k].rgb[2] = 0.01
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].wa = 0
			lp[k].wb = 90.0
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			lp[k].P[0].X = tree[i].x - (tree[i].W1 * 0.5)
			lp[k].P[0].Y = tree[i].y - (tree[i].W1 * 0.5)
			lp[k].P[0].Z = tree[i].z
			lp[k].P[1].X = lp[k].P[0].X
			lp[k].P[1].Y = lp[k].P[0].Y
			lp[k].P[1].Z = tree[i].z + tree[i].H1
			lp[k].P[2].X = tree[i].x + (tree[i].W1 * 0.5)
			lp[k].P[2].Y = lp[k].P[0].Y
			lp[k].P[2].Z = tree[i].z + tree[i].H1
			lp[k].P[3].X = lp[k].P[2].X
			lp[k].P[3].Y = lp[k].P[0].Y
			lp[k].P[3].Z = tree[i].z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1

			/*----2----*/
			name = tree[i].treename + "-m2"
			lp[k].opname = name
			lp[k].wa = -90.0
			lp[k].wb = 90.0
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.4
			lp[k].rgb[1] = 0.3
			lp[k].rgb[2] = 0.01
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x + (tree[i].W1 * 0.5)
			lp[k].P[2].Y = tree[i].y + (tree[i].W1 * 0.5)
			lp[k].P[2].Z = tree[i].z + tree[i].H1
			lp[k].P[3].X = lp[k].P[2].X
			lp[k].P[3].Y = lp[k].P[2].Y
			lp[k].P[3].Z = tree[i].z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1

			/*----3----*/
			name = tree[i].treename + "-m3"
			lp[k].opname = name

			lp[k].wa = 180.0
			lp[k].wb = 90.0
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.4
			lp[k].rgb[1] = 0.3
			lp[k].rgb[2] = 0.01
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W1 * 0.5)
			lp[k].P[2].Y = tree[i].y + (tree[i].W1 * 0.5)
			lp[k].P[2].Z = tree[i].z + tree[i].H1
			lp[k].P[3].X = lp[k].P[2].X
			lp[k].P[3].Y = lp[k].P[2].Y
			lp[k].P[3].Z = tree[i].z

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1

			/*----4----*/
			name = tree[i].treename + "-m4"
			lp[k].opname = name

			lp[k].wa = 90.0
			lp[k].wb = 90.0
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.4
			lp[k].rgb[1] = 0.3
			lp[k].rgb[2] = 0.01
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2] = lp[k-3].P[1]
			lp[k].P[3] = lp[k-3].P[0]

			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)
			//HOUSEN2(&lp[k].P[0],&lp[k].P[1],&lp[k].P[2],&lp[k].e)

			k = k + 1

			/*----5----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = -22.5
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0].X = tree[i].x
			lp[k].P[0].Y = tree[i].y - (tree[i].W2 * 0.5)
			lp[k].P[0].Z = tree[i].z + tree[i].H1
			lp[k].P[1].X = tree[i].x
			lp[k].P[1].Y = tree[i].y - (tree[i].W3 * 0.5)
			lp[k].P[1].Z = tree[i].z + tree[i].H1 + tree[i].H2
			lp[k].P[2].X = tree[i].x + (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y - (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x + (tree[i].W2*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y - (tree[i].W2*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = tree[i].z + tree[i].H1

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----6----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x + (tree[i].W3 * 0.5)
			lp[k].P[2].Y = tree[i].y
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x + (tree[i].W2 * 0.5)
			lp[k].P[3].Y = tree[i].y
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----7----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x + (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y + (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x + (tree[i].W2*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y + (tree[i].W2*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----8----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x
			lp[k].P[2].Y = tree[i].y + (tree[i].W3 * 0.5)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x
			lp[k].P[3].Y = tree[i].y + (tree[i].W2 * 0.5)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----9----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = 360 + (lp[k-1].wa - 45)
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y + (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x - (tree[i].W2*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y + (tree[i].W2*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----10----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W3 * 0.5)
			lp[k].P[2].Y = tree[i].y
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x - (tree[i].W2 * 0.5)
			lp[k].P[3].Y = tree[i].y
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----11----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y - (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x - (tree[i].W2*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y - (tree[i].W2*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----12----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x
			lp[k].P[2].Y = tree[i].y - (tree[i].W3 * 0.5)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x
			lp[k].P[3].Y = tree[i].y - (tree[i].W2 * 0.5)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----13----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = -22.5
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0].X = tree[i].x
			lp[k].P[0].Y = tree[i].y - (tree[i].W3 * 0.5)
			lp[k].P[0].Z = tree[i].z + tree[i].H1 + tree[i].H2
			lp[k].P[1].X = tree[i].x
			lp[k].P[1].Y = tree[i].y - (tree[i].W4 * 0.5)
			lp[k].P[1].Z = tree[i].z + tree[i].H1 + tree[i].H2 + tree[i].H3
			lp[k].P[2].X = tree[i].x + (tree[i].W4*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y - (tree[i].W4*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x + (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y - (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----14----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x + (tree[i].W4 * 0.5)
			lp[k].P[2].Y = tree[i].y
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x + (tree[i].W3 * 0.5)
			lp[k].P[3].Y = tree[i].y
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----15----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x + (tree[i].W4*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y + (tree[i].W4*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x + (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y + (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----16----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x
			lp[k].P[2].Y = tree[i].y + (tree[i].W4 * 0.5)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x
			lp[k].P[3].Y = tree[i].y + (tree[i].W3 * 0.5)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*---17---*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = 360 + (lp[k-1].wa - 45)
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W4*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y + (tree[i].W4*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x - (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y + (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*----18----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W4 * 0.5)
			lp[k].P[2].Y = tree[i].y
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x - (tree[i].W3 * 0.5)
			lp[k].P[3].Y = tree[i].y
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*---19----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x - (tree[i].W4*0.5)*math.Cos(45*M_rad)
			lp[k].P[2].Y = tree[i].y - (tree[i].W4*0.5)*math.Sin(45*M_rad)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x - (tree[i].W3*0.5)*math.Cos(45*M_rad)
			lp[k].P[3].Y = tree[i].y - (tree[i].W3*0.5)*math.Sin(45*M_rad)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1

			/*---20----*/
			name = tree[i].treename
			lp[k].opname = name

			lp[k].wa = lp[k-1].wa - 45
			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			lp[k].rgb[0] = 0.0
			lp[k].rgb[1] = 1
			lp[k].rgb[2] = 0.0
			lp[k].polyd = 4
			lp[k].P = make([]XYZ, 4)

			lp[k].P[0] = lp[k-1].P[3]
			lp[k].P[1] = lp[k-1].P[2]
			lp[k].P[2].X = tree[i].x
			lp[k].P[2].Y = tree[i].y - (tree[i].W4 * 0.5)
			lp[k].P[2].Z = lp[k].P[1].Z
			lp[k].P[3].X = tree[i].x
			lp[k].P[3].Y = tree[i].y - (tree[i].W3 * 0.5)
			lp[k].P[3].Z = lp[k].P[0].Z

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			k = k + 1
		} else {
			fmt.Printf("error--**COORDNT-lp TREE \n")
			os.Exit(1)
		}
	}

	/*-------多角形の障害物の直接入力----------------------*/
	for i = 0; i < polyn; i++ {

		if poly[i].polyknd == "OBS" {
			lp[k].opname = poly[i].polyname

			lp[k].polyd = poly[i].polyd
			lp[k].P = make([]XYZ, lp[k].polyd)

			for j = 0; j < lp[k].polyd; j++ {
				lp[k].P[j] = poly[i].P[j]
			}

			HOUSEN2(&lp[k].P[0], &lp[k].P[1], &lp[k].P[2], &lp[k].e)
			lp[k].wb = math.Acos(lp[k].e.Z) * (180 / math.Pi)
			lp[k].wa = math.Asin(lp[k].e.X/math.Sin(lp[k].wb*M_rad)) * (180 / math.Pi)
			/*--法線ベクトルの算出　HOUSEN2関数を使うと、向きが逆になるので、変更 091128 higuchi --*/
			lp[k].e.Z = math.Cos(lp[k].wb * M_rad)
			lp[k].e.Y = -math.Sin(lp[k].wb*M_rad) * math.Cos(lp[k].wa*M_rad)
			lp[k].e.X = -math.Sin(lp[k].wb*M_rad) * math.Sin(lp[k].wa*M_rad)
			CAT(&lp[k].e.X, &lp[k].e.Y, &lp[k].e.Z)

			lp[k].ref = 0.0
			lp[k].wd = 0
			lp[k].sbflg = 0

			lp[k].rgb[0] = poly[i].rgb[0]
			lp[k].rgb[1] = poly[i].rgb[1]
			lp[k].rgb[2] = poly[i].rgb[2]

			for p = 0; p < 366; p++ {
				lp[k].shad[p] = 1
			}

			k = k + 1
		}
	}
	*lpn = k
}

/*----------------------

受照面OPに変換する関数               1999.5.22
------------------------*/
func OP_COORDNT(opn *int, bdpn int, BDP []BBDP, op []P_MENN, polyn int, poly []POLYGN) {
	var l, j, i, v, sumop, k int
	var cWa, sWa, cWb, sWb, a, b, c float64
	var ax, by, cz float64

	sumop = 0

	const M_rad = math.Pi / 180.0

	for l = 0; l < bdpn; l++ {
		k = 0
		cWa = 0.0
		sWa = 0.0
		cWb = 0.0
		sWb = 0.0

		cWa = math.Cos(-BDP[l].Wa * M_rad)
		sWa = math.Sin(-BDP[l].Wa * M_rad)
		cWb = math.Cos(BDP[l].Wb * M_rad)
		sWb = math.Sin(BDP[l].Wb * M_rad)

		sumop = sumop + BDP[l].sumRMP
		v = sumop - BDP[l].sumRMP
		for j = v; j < sumop; j++ {
			a = 0.0
			b = 0.0
			c = 0.0

			op[j].wd = BDP[l].RMP[k].sumWD
			op[j].wlflg = 0

			op[j].opname = BDP[l].RMP[k].rmpname

			op[j].rgb[0] = BDP[l].RMP[k].rgb[0]
			op[j].rgb[1] = BDP[l].RMP[k].rgb[1]
			op[j].rgb[2] = BDP[l].RMP[k].rgb[2]
			op[j].polyd = 4
			op[j].P = make([]XYZ, 4)

			a = BDP[l].x0 + BDP[l].RMP[k].xb0*cWa - BDP[l].RMP[k].yb0*cWb*sWa
			b = BDP[l].y0 + BDP[l].RMP[k].xb0*sWa + BDP[l].RMP[k].yb0*cWb*cWa
			c = BDP[l].z0 + BDP[l].RMP[k].yb0*sWb

			op[j].P[0].X = a
			op[j].P[0].Y = b
			op[j].P[0].Z = c

			op[j].P[1].X = a - BDP[l].RMP[k].Rh*cWb*sWa
			op[j].P[1].Y = b + BDP[l].RMP[k].Rh*cWb*cWa
			op[j].P[1].Z = c + BDP[l].RMP[k].Rh*sWb

			op[j].P[2].X = BDP[l].RMP[k].Rw*cWa + op[j].P[1].X
			op[j].P[2].Y = BDP[l].RMP[k].Rw*sWa + op[j].P[1].Y
			op[j].P[2].Z = op[j].P[1].Z

			op[j].P[3].X = a + BDP[l].RMP[k].Rw*cWa
			op[j].P[3].Y = b + BDP[l].RMP[k].Rw*sWa
			op[j].P[3].Z = c

			op[j].grpx = BDP[l].RMP[k].grpx
			op[j].wa = BDP[l].Wa
			op[j].wb = BDP[l].Wb
			op[j].ref = BDP[l].RMP[k].ref
			op[j].e.Z = math.Cos(op[j].wb * M_rad)
			op[j].e.Y = -math.Sin(op[j].wb*M_rad) * math.Cos(op[j].wa*M_rad)
			op[j].e.X = -math.Sin(op[j].wb*M_rad) * math.Sin(op[j].wa*M_rad)
			CAT(&op[j].e.X, &op[j].e.Y, &op[j].e.Z)

			op[j].Nopw = BDP[l].RMP[k].sumWD
			// opw構造体のメモリ確保
			if op[j].Nopw > 0 {
				op[j].opw = make([]WD_MENN, op[j].Nopw)
				for m := 0; m < op[j].Nopw; m++ {
					opw := op[j].opw[m]

					opw.opwname = ""
					opw.P = nil
					opw.polyd = 0
					opw.ref = 0.0
					opw.grpx = 0.0
					opw.sumw = 0.0
					matinit(opw.rgb[:], 3)
				}
			}
			if BDP[l].RMP[k].sumWD != 0 {
				for i = 0; i < BDP[l].RMP[k].sumWD; i++ {
					ax = 0.0
					by = 0.0
					cz = 0.0

					op[j].opw[i].opwname = BDP[l].RMP[k].WD[i].winname

					ax = a + BDP[l].RMP[k].WD[i].xr*cWa - BDP[l].RMP[k].WD[i].yr*cWb*sWa
					by = b + BDP[l].RMP[k].WD[i].xr*sWa + BDP[l].RMP[k].WD[i].yr*cWb*cWa
					cz = c + BDP[l].RMP[k].WD[i].yr*sWb

					op[j].opw[i].grpx = BDP[l].RMP[k].WD[i].grpx
					op[j].opw[i].ref = BDP[l].RMP[k].WD[i].ref

					op[j].opw[i].rgb[0] = BDP[l].RMP[k].WD[i].rgb[0]
					op[j].opw[i].rgb[1] = BDP[l].RMP[k].WD[i].rgb[1]
					op[j].opw[i].rgb[2] = BDP[l].RMP[k].WD[i].rgb[2]
					op[j].opw[i].polyd = 4
					op[j].opw[i].P = make([]XYZ, 4)

					op[j].opw[i].P[0].X = ax
					op[j].opw[i].P[0].Y = by
					op[j].opw[i].P[0].Z = cz

					op[j].opw[i].P[1].X = ax - BDP[l].RMP[k].WD[i].Wh*cWb*sWa
					op[j].opw[i].P[1].Y = by + BDP[l].RMP[k].WD[i].Wh*cWb*cWa
					op[j].opw[i].P[1].Z = cz + BDP[l].RMP[k].WD[i].Wh*sWb
					op[j].opw[i].P[2].X = ax + BDP[l].RMP[k].WD[i].Ww*cWa - BDP[l].RMP[k].WD[i].Wh*cWb*sWa
					op[j].opw[i].P[2].Y = by + BDP[l].RMP[k].WD[i].Ww*sWa + BDP[l].RMP[k].WD[i].Wh*cWb*cWa
					op[j].opw[i].P[2].Z = cz + BDP[l].RMP[k].WD[i].Wh*sWb
					op[j].opw[i].P[3].X = ax + BDP[l].RMP[k].WD[i].Ww*cWa
					op[j].opw[i].P[3].Y = by + BDP[l].RMP[k].WD[i].Ww*sWa
					op[j].opw[i].P[3].Z = cz
				}
			}
			k = k + 1
		}
	}

	for i = 0; i < polyn; i++ {
		if poly[i].polyknd == "RMP" {

			op[sumop].opname = poly[i].polyname
			op[sumop].ref = poly[i].ref
			op[sumop].wd = 0
			op[sumop].grpx = poly[i].grpx

			op[sumop].rgb[0] = poly[i].rgb[0]
			op[sumop].rgb[1] = poly[i].rgb[1]
			op[sumop].rgb[2] = poly[i].rgb[2]
			op[sumop].polyd = poly[i].polyd
			op[sumop].P = make([]XYZ, op[sumop].polyd)
			for j = 0; j < op[sumop].polyd; j++ {
				op[sumop].P[j] = poly[i].P[j]
			}

			HOUSEN2(&op[sumop].P[0], &op[sumop].P[1], &op[sumop].P[2], &op[sumop].e)
			op[sumop].wb = math.Acos(op[sumop].e.Z) * (180 / math.Pi)
			op[sumop].wa = math.Asin(op[sumop].e.X/math.Sin(op[sumop].wb*M_rad)) * (180 / math.Pi)

			/*--ポリゴンの法線ベクトルが逆で出てしまうのを修正 091128 higuchi--*/
			op[sumop].e.Z = math.Cos(op[sumop].wb * M_rad)
			op[sumop].e.Y = -math.Sin(op[sumop].wb*M_rad) * math.Cos(op[sumop].wa*M_rad)
			op[sumop].e.X = -math.Sin(op[sumop].wb*M_rad) * math.Sin(op[sumop].wa*M_rad)

			sumop = sumop + 1
		}

	}

	*opn = sumop
}