User:Jorge Stolfi/make-coord-system-figure
Jump to navigation
Jump to search
#! /usr/bin/python -t # _*_ coding: iso-8859-1 _*_ # Last edited on 2009-05-03 18:12:55 by stolfilocal PROG_NAME = "make-coord-system-figure" PROG_DESC = "Generates an SVG illustration for the Wikipedia articles on coord systems" PROG_VERS = "1.0" import sys import re import os import copy import math from math import sqrt,sin,cos sys.path[1:0] = [ sys.path[0] + '/../lib', os.path.expandvars('${STOLFIHOME}/lib'), '.' ] import argparser; from argparser import ArgParser import rn import rmxn import hrn import perspective from decimal import * from datetime import date PROG_COPYRIGHT = "Copyright © 2009-05-02 by the State University of Campinas (UNICAMP)" PROG_HELP = \ PROG_NAME+ " \\\n" \ " -back {BOOL} \\\n" \ " -frame {BOOL} \\\n" \ " -rho {BOOL} \\\n" \ " -system { \"SE\" | \"SZ\" | \"CY\" | \"CA\" } \\\n" \ +argparser.help_info_HELP+ " \\\n" \ " > {FIGURE}.svg" PROG_INFO = \ "NAME\n" \ " " +PROG_NAME+ " - " +PROG_DESC+ ".\n" \ "\n" \ "SYNOPSIS\n" \ " " +PROG_HELP+ "\n" \ "\n" \ "DESCRIPTION\n" \ " Writes an SVG illustration for the Wikipedia articles on coord systems.\n" \ "\n" \ "OPTIONS\n" \ " -back {BOOL} \n" \ " If {BOOL} is 1, paints the background with a nonwhite color. If {BOOL} is 0," \ " leaves it transparent. This option is meant to debug the image size.\n" \ "\n" \ " -frame {BOOL} \n" \ " If {BOOL} is 1, draws some framing lines. This option is meant" \ " to debug the plot dimensions.\n" \ "\n" \ " -rho {BOOL} \n" \ " If {BOOL} is 1, uses \"rho\" for radius, else uses \"r\".\n" \ "\n" \ " -system { \"SE\" | \"SZ\" | \"CY\" | \"CA\" } \n" \ " Coordinate system to illustrate.\n" \ " \"SE\" = Spherical (elevation angle).\n" \ " \"SZ\" = Spherical (zenith angle).\n" \ " \"CY\" = Cylindrical.\n" \ " \"CA\" = Cartesian.\n" \ "\n" \ "DOCUMENTATION OPTIONS\n" \ +argparser.help_info_INFO+ "\n" \ "\n" \ "SEE ALSO\n" \ " cat(1).\n" \ "\n" \ "AUTHOR\n" \ " Created 2009-04-04 by Jorge Stolfi, IC-UNICAMP.\n" \ "\n" \ "MODIFICATION HISTORY\n" \ " 2009-04-04 by J. Stolfi, IC-UNICAMP: created.\n" \ "\n" \ "WARRANTY\n" \ " " +argparser.help_info_NO_WARRANTY+ "\n" \ "\n" \ "RIGHTS\n" \ " " +PROG_COPYRIGHT+ ".\n" \ "\n" \ " " +argparser.help_info_STANDARD_RIGHTS # COMMAND ARGUMENT PARSING pp = ArgParser(sys.argv, sys.stderr, PROG_HELP, PROG_INFO) class Options : back = None; system = None; err = None; def arg_error(msg): "Prints the error message {msg} about the command line arguments, and aborts." sys.stderr.write("%s\n" % msg); sys.stderr.write("usage: %s\n" % PROG_HELP); sys.exit(1) def parse_args(pp) : "Parses command line arguments.\n" \ "\n" \ " Expects an {ArgParser} instance {pp} containing the arguments," \ " still unparsed. Returns an {Options} instance {op}, where" \ " {op.err} is an error message, if any (a string) or {None}." op = Options(); # Being optimistic: op.err = None pp.get_keyword("-back") op.back = pp.get_next_int(0, 1) pp.get_keyword("-frame") op.frame = pp.get_next_int(0, 1) pp.get_keyword("-rho") op.rho = pp.get_next_int(0, 1) pp.get_keyword("-system") op.system = pp.get_next() return op #---------------------------------------------------------------------- class Dimensions : "Plot dimensions and perspective matrix." def __init__(dim, op) : dim.map = None; # World to Image perspective map. # Coordinates of point: dim.xc_val = None; # X coordinate (for all systems). dim.yc_val = None; # Y coordinate (for all systems). dim.zc_val = None; # Z coordinate (for all systems). dim.rc_val = None; # Radial coordinate (for SE, SZ, CY). dim.ac_val = None; # Azimuth angle (for SE, SZ, CY). dim.ec_val = None; # Elevation angle (for SE, SZ). dim.hd_len = 0.08; # Length of arrow heads (W). dim.ax_len = 1.20; # Length of coord axes (W). dim.font_wy = 30; # Font size in pixels. dim.pt_rad = 0.02; # Point size in World units. dim.ax_unit = 0.2; # Unit for axis ticks. dim.fig_wx = 620; # Total figure width. dim.fig_wy = 600; # Total figure height. dim.ct_color = '200,0,100'; # Color for significant coord traces. dim.map = None; # World to Image perspective map. dim.scale = 1.0; # Final scale factor (can be chaged without changing any other dim). # Coordinates of point: if (op.system == 'CA') : dim.xc_val = 0.4; dim.yc_val = 0.6; dim.zc_val = 0.8; elif (op.system == 'CY') : dim.rc_val = 0.8; dim.ac_val = math.radians(130); dim.zc_val = 0.8; # Compute X,Y from azimuth and radius: dim.xc_val = dim.rc_val*cos(dim.ac_val); dim.yc_val = dim.rc_val*sin(dim.ac_val); elif (op.system == 'SE'): dim.rc_val = 0.8; dim.ac_val = math.radians(130); dim.ec_val = math.radians(50); # Compute X,Y,Z from azimuth, elevation, and radius: dim.xc_val = dim.rc_val*cos(dim.ec_val)*cos(dim.ac_val); dim.yc_val = dim.rc_val*cos(dim.ec_val)*sin(dim.ac_val); dim.zc_val = dim.rc_val*sin(dim.ec_val); elif (op.system == 'SZ'): dim.rc_val = 0.8; dim.ac_val = math.radians(130); dim.ec_val = math.radians(70); # Compute X,Y,Z from azimuth, elevation, and radius: dim.xc_val = dim.rc_val*sin(dim.ec_val)*cos(dim.ac_val); dim.yc_val = dim.rc_val*sin(dim.ec_val)*sin(dim.ac_val); dim.zc_val = dim.rc_val*cos(dim.ec_val); else : assert False, "invalid coord system"; # Perspective map: att = [0,0,0]; obs = [5,3,3]; upd = [0,0,1]; mag = [+220,-220,+220]; ctr = [+320,+300,0000]; dim.map = perspective.camera_matrix(att,obs,upd,mag,ctr); sys.stderr.write("dim.map = %s\n" % dim.map); #---------------------------------------------------------------------- #---------------------------------------------------------------------- def output_figure(op) : "Writes the figure to {stdout}." "\n" \ " Expects an {Options} instance {op}." # Computes the sizes of things and the perspective map: dim = Dimensions(op) output_figure_preamble(op,dim) sys.stdout.write('\n') output_figure_obj_defs(op,dim) sys.stdout.write('\n') output_figure_body(op,dim) sys.stdout.write('\n') output_figure_postamble(op,dim) sys.stderr.write("done.\n") #---------------------------------------------------------------------- def output_figure_preamble(op,dim) : "Writes the SVG preamble to {stdout}." sys.stdout.write( \ '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n' \ '<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">\n' \ '<!-- Created on '+ date.isoformat(date.today()) +' by Jorge Stolfi with the script make-coord-system-figure -->\n' \ '<!-- This file is declared PUBLIC DOMAIN by its creator -->\n' \ '\n' \ '<svg\n' \ ' id="fig"\n' \ ' xmlns="http://www.w3.org/2000/svg"\n' \ ' xmlns:xlink="http://www.w3.org/1999/xlink"\n' \ '\n' \ ' fill="none"\n' \ ' fill-opacity="1"\n' \ ' fill-rule="evenodd"\n' \ '\n' \ ' stroke-linecap="round"\n' \ ' stroke-linejoin="round"\n' \ ' stroke-dasharray="none"\n' \ ' stroke-opacity="1"\n' \ ' stroke-width="1.5px"\n' \ '\n' \ ' font-style="normal"\n' \ ' font-weight="normal"\n' \ ' font-size="'+ dts(dim.font_wy) +'px"\n' \ ' font-family="Bitstream Vera"\n' \ '\n' \ ' width="'+ dts(dim.scale*dim.fig_wx) +'"\n' \ ' height="'+ dts(dim.scale*dim.fig_wy) +'"\n' \ '>\n' \ ) sys.stdout.write('\n') if (op.back) : sys.stdout.write( \ ' <!-- Rectangle to test the figure size -->\n' \ ' <rect x="000" y="000" width="'+ dts(dim.scale*dim.fig_wx) +'" height="'+ dts(dim.scale*dim.fig_wy) +'"' \ ' stroke="#000000" stroke-width="2px" fill="#ffcc55"' \ ' />\n' \ ) sys.stdout.write('\n') # Scale everything and position the diagram: sys.stdout.write( \ ' <g\n' \ ' transform="scale('+ dts(dim.scale)+ ')"\n' \ ' >\n' \ ) #---------------------------------------------------------------------- def output_figure_obj_defs(op,dim) : "Writes the main object definitions to {stdout}." sys.stdout.write( \ ' <defs>\n' \ ' </defs>\n' \ ) #---------------------------------------------------------------------- def output_figure_body(op,dim) : "Writes the body of the figure to {stdout}." # Plot the reference plane(s): if (op.system == 'CA') : output_system_planes_CA(op,dim); else: output_system_planes_CY_SE_SZ(op,dim); pos_W = [dim.xc_val,dim.yc_val,dim.zc_val]; if (op.system == 'SE') : output_system_SE(op,dim,pos_W); elif (op.system == 'CA') : output_system_CA(op,dim,pos_W); elif (op.system == 'CY') : output_system_CY(op,dim,pos_W); elif (op.system == 'SZ') : output_system_SZ(op,dim,pos_W); else : assert False, "invalid {op.system}"; def output_system_planes_CA(op,dim): # Key World points: ooo_W = [00,00,00]; # Origin. poo_W = [+1,00,00]; # +X axis dir. moo_W = [-1,00,00]; # -X axis dir. opo_W = [00,+1,00]; # +Y axis dir. omo_W = [00,-1,00]; # -Y axis dir. oop_W = [00,00,+1]; # +Z axis dir. oom_W = [00,00,-1]; # -Z axis dir. # Key Image points: ooo = img_point(ooo_W,dim.map); # Origin. moo = img_point(moo_W,dim.map); # Point -1 on X axis poo = img_point(poo_W,dim.map); # Point +1 on X axis omo = img_point(omo_W,dim.map); # Point -1 on U axis opo = img_point(opo_W,dim.map); # Point +1 on U axis oom = img_point(oom_W,dim.map); # Point -1 on Z axis oop = img_point(oop_W,dim.map); # Point +1 on Z axis pmo = img_point(rn.add(poo_W,omo_W),dim.map); # Corner A of XY square. ppo = img_point(rn.add(poo_W,opo_W),dim.map); # Corner B of XY square. mpo = img_point(rn.add(moo_W,opo_W),dim.map); # Corner C of XY square. mmo = img_point(rn.add(moo_W,omo_W),dim.map); # Corner D of XY square. opm = img_point(rn.add(opo_W,oom_W),dim.map); # Corner A of UZ square. opp = img_point(rn.add(opo_W,oop_W),dim.map); # Corner B of UZ square. omp = img_point(rn.add(omo_W,oop_W),dim.map); # Corner C of UZ square. omm = img_point(rn.add(omo_W,oom_W),dim.map); # Corner D of UZ square. mop = img_point(rn.add(moo_W,oop_W),dim.map); # Corner A of ZX square. pop = img_point(rn.add(poo_W,oop_W),dim.map); # Corner B of ZX square. pom = img_point(rn.add(poo_W,oom_W),dim.map); # Corner C of ZX square. mom = img_point(rn.add(moo_W,oom_W),dim.map); # Corner D of ZX square. sys.stdout.write( \ ' <!-- The reference planes: -->\n' \ ' <g stroke="rgb(185, 205, 170)" fill="rgb(215, 235, 200)" fill-opacity="0.5">\n' \ ' <polygon points="'+ pfm(moo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(mom) +'"/>\n' \ ' <polygon points="'+ pfm(omo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(omm) +'"/>\n' \ ' <polygon points="'+ pfm(poo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(pom) +'"/>\n' \ ' <polygon points="'+ pfm(opo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(opm) +'"/>\n' \ ' \n' \ ' <polygon points="'+ pfm(pmo) +' '+ pfm(ppo) +' '+ pfm(mpo) +' '+ pfm(mmo) +'"/>\n' \ ' \n' \ ' <polygon points="'+ pfm(moo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(mop) +'"/>\n' \ ' <polygon points="'+ pfm(omo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(omp) +'"/>\n' \ ' <polygon points="'+ pfm(poo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(pop) +'"/>\n' \ ' <polygon points="'+ pfm(opo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(opp) +'"/>\n' \ ' </g>\n' \ ); if (op.frame) : # Auxiliary lines for debugging: sys.stdout.write( \ ' <polygon points="'+ pfm(pmo) +' '+ pfm(ppo) +' '+ pfm(mpo) +' '+ pfm(mmo) +'" stroke="rgb(177,0,0)" fill="none" />\n' \ ' <path d="M '+ pfm(pmo) +' L '+ pfm(mpo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ' <path d="M '+ pfm(ppo) +' L '+ pfm(mmo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ' <path d="M '+ pfm(moo) +' L '+ pfm(poo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ' <path d="M '+ pfm(omo) +' L '+ pfm(opo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ); sys.stdout.write('\n'); #---------------------------------------------------------------------- def output_system_planes_CY_SE_SZ(op,dim): # Key World points: ooo_W = [00,00,00]; # Origin. poo_W = [+1,00,00]; # +X axis dir. moo_W = [-1,00,00]; # -X axis dir. opo_W = vec_from_ang_rad([+1,00,00],[00,+1,00],dim.ac_val,1.0); # +U axis dir. omo_W = vec_from_ang_rad([-1,00,00],[00,-1,00],dim.ac_val,1.0); # -U axis dir. oop_W = [00,00,+1]; # +Z axis dir. oom_W = [00,00,-1]; # -Z axis dir. # Key Image points: ooo = img_point(ooo_W,dim.map); # Origin. moo = img_point(moo_W,dim.map); # Point -1 on X axis poo = img_point(poo_W,dim.map); # Point +1 on X axis omo = img_point(omo_W,dim.map); # Point -1 on U axis opo = img_point(opo_W,dim.map); # Point +1 on U axis oom = img_point(oom_W,dim.map); # Point -1 on Z axis oop = img_point(oop_W,dim.map); # Point +1 on Z axis opm = img_point(rn.add(opo_W,oom_W),dim.map); # Corner A of UZ square. opp = img_point(rn.add(opo_W,oop_W),dim.map); # Corner B of UZ square. omp = img_point(rn.add(omo_W,oop_W),dim.map); # Corner C of UZ square. omm = img_point(rn.add(omo_W,oom_W),dim.map); # Corner D of UZ square. mop = img_point(rn.add(moo_W,oop_W),dim.map); # Corner A of ZX square. pop = img_point(rn.add(poo_W,oop_W),dim.map); # Corner B of ZX square. pom = img_point(rn.add(poo_W,oom_W),dim.map); # Corner C of ZX square. mom = img_point(rn.add(moo_W,oom_W),dim.map); # Corner D of ZX square. # !!! BUG !!! should compute the ellipse from the map !!! sys.stdout.write( \ ' <!-- The reference planes: -->\n' \ ' <g stroke="rgb(185, 205, 170)" fill="rgb(215, 235, 200)" fill-opacity="0.5">\n' \ ' <polygon points="'+ pfm(moo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(mom) +'"/>\n' \ ' <polygon points="'+ pfm(omo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(omm) +'"/>\n' \ ' <polygon points="'+ pfm(poo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(pom) +'"/>\n' \ ' <polygon points="'+ pfm(opo) +' '+ pfm(ooo) +' '+ pfm(oom) +' '+ pfm(opm) +'"/>\n' \ ' \n' \ ' <ellipse cx="0" cy="14" rx="222" ry="102" transform="translate('+ pfm(ooo) +')rotate(0)"/>\n' \ ' \n' \ ' <polygon points="'+ pfm(moo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(mop) +'"/>\n' \ ' <polygon points="'+ pfm(omo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(omp) +'"/>\n' \ ' <polygon points="'+ pfm(poo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(pop) +'"/>\n' \ ' <polygon points="'+ pfm(opo) +' '+ pfm(ooo) +' '+ pfm(oop) +' '+ pfm(opp) +'"/>\n' \ ' </g>\n' \ ); sys.stdout.write('\n'); if (op.frame) : # Auxiliary lines for debugging: pmo = img_point([+1,-1,00],dim.map); # Corner A of XY square. ppo = img_point([+1,+1,00],dim.map); # Corner B of XY square. mpo = img_point([-1,+1,00],dim.map); # Corner C of XY square. mmo = img_point([-1,-1,00],dim.map); # Corner D of XY square. sys.stdout.write( \ ' <polygon points="'+ pfm(pmo) +' '+ pfm(ppo) +' '+ pfm(mpo) +' '+ pfm(mmo) +'" stroke="rgb(177,0,0)" fill="none" />\n' \ ' <path d="M '+ pfm(pmo) +' L '+ pfm(mpo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ' <path d="M '+ pfm(ppo) +' L '+ pfm(mmo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ' <path d="M '+ pfm(moo) +' L '+ pfm(poo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ' <path d="M '+ pfm(omo) +' L '+ pfm(opo) +'" stroke="rgb(137, 137, 137)"/>\n' \ ); sys.stdout.write('\n'); #---------------------------------------------------------------------- def output_system_CA(op,dim,pos_W) : # Relevant points: ooo_W = [00,00,00]; px_W = [dim.xc_val,00,00]; py_W = [00,dim.yc_val,00]; pz_W = [00,00,dim.zc_val]; pxy_W = [dim.xc_val,dim.yc_val,00]; pyz_W = [00,dim.yc_val,dim.zc_val]; pxz_W = [dim.xc_val,00,dim.zc_val]; output_coord_axis(op,dim,[+1,00,00],[-10,+5],'X'); output_coord_axis(op,dim,[00,+1,00],[+5,-5],'Y'); output_coord_axis(op,dim,[00,00,+1],[-7,+5],'Z'); sys.stdout.write('\n'); ooo = img_point(ooo_W,dim.map) output_label(op,dim,ooo,[-5,-5],make_italic_style('O')); # Coordinate traces: output_straight_trace(op,dim,ooo_W,px_W,True,'x'); # output_straight_trace(op,dim,ooo_W,py_W,False,''); # output_straight_trace(op,dim,ooo_W,pz_W,False,''); output_straight_trace(op,dim,px_W,pxy_W,True,'y'); output_straight_trace(op,dim,px_W,pxz_W,True,''); output_straight_trace(op,dim,py_W,pxy_W,True,''); output_straight_trace(op,dim,py_W,pyz_W,True,''); output_straight_trace(op,dim,pz_W,pxz_W,True,''); output_straight_trace(op,dim,pz_W,pyz_W,True,''); output_straight_trace(op,dim,pyz_W,pos_W,True,''); output_straight_trace(op,dim,pxz_W,pos_W,True,''); output_straight_trace(op,dim,pxy_W,pos_W,True,'z'); # The point: output_point(op,dim,pos_W,dim.pt_rad,'x','y','z'); #---------------------------------------------------------------------- def output_system_CY(op,dim,pos_W) : output_coord_axis(op,dim,[+1,00,00],[-5,+5],'A'); output_coord_axis(op,dim,[00,00,+1],[-5,-5],'L'); sys.stdout.write('\n'); if (op.rho) : r_str = 'ρ'; else : r_str = 'r'; # Relevant points: ooo_W = [00,00,00]; Pr_W = [dim.rc_val,00,00]; Pz_W = [00,00,dim.zc_val]; Pra_W = vec_from_ang_rad([+1,00,00],[00,+1,00],dim.ac_val,dim.rc_val); Prz_W = [dim.rc_val,00,dim.zc_val]; ooo = img_point(ooo_W,dim.map) output_label(op,dim,ooo,[-5,-5],make_italic_style('O')); # Radial traces: output_straight_trace(op,dim,ooo_W,Pr_W,True,r_str); output_straight_trace(op,dim,ooo_W,Pra_W,True,''); output_straight_trace(op,dim,Pz_W,Prz_W,True,''); output_straight_trace(op,dim,Pz_W,pos_W,True,''); # Azimuth traces: output_circular_trace(op,dim,ooo_W,[+1,00,00],[00,+1,00],0,dim.ac_val,dim.rc_val,'φ'); output_circular_trace(op,dim,Pz_W,[+1,00,00],[00,+1,00],0,dim.ac_val,dim.rc_val,''); # Longitudinal traces: # output_straight_trace(op,dim,ooo_W,Pz_W,False,''); output_straight_trace(op,dim,Pr_W,Prz_W,True,''); output_straight_trace(op,dim,Pra_W,pos_W,True,'z'); # The point: output_point(op,dim,pos_W,dim.pt_rad,r_str,'φ','z'); #---------------------------------------------------------------------- def output_system_SE(op,dim,pos_W) : output_coord_axis(op,dim,[+1,00,00],[-5,+5],'A'); sys.stdout.write('\n'); if (op.rho) : r_str = 'ρ'; else : r_str = 'r'; # Relevant points: ooo_W = [00,00,00]; Pr_W = [dim.rc_val,00,00]; Pz_W = [00,00,dim.zc_val]; Pra_W = vec_from_ang_rad([+1,00,00],[00,+1,00],dim.ac_val,dim.rc_val); Pre_W = vec_from_ang_rad([+1,00,00],[00,00,+1],dim.ec_val,dim.rc_val); # Radial traces: output_straight_trace(op,dim,ooo_W,Pr_W,True,r_str); output_straight_trace(op,dim,ooo_W,Pra_W,True,''); output_straight_trace(op,dim,ooo_W,Pre_W,True,''); output_straight_trace(op,dim,ooo_W,pos_W,True,''); output_straight_trace(op,dim,Pz_W,Pre_W,True,''); output_straight_trace(op,dim,Pz_W,pos_W,True,''); output_straight_trace(op,dim,ooo_W,Pz_W,True,''); # Elevation traces: xy_W = vec_from_ang_rad([+1,00,00],[00,+1,00],dim.ac_val,1.0); output_circular_trace(op,dim,ooo_W,[+1,00,00],[00,00,+1],0,dim.ec_val,dim.rc_val,'θ'); output_circular_trace(op,dim,ooo_W,xy_W,[00,00,+1],0,dim.ec_val,dim.rc_val,''); # Azimuth traces: tp_rad = dim.rc_val*cos(dim.ec_val); output_circular_trace(op,dim,ooo_W,[+1,00,00],[00,+1,00],0,dim.ac_val,dim.rc_val,''); output_circular_trace(op,dim,Pz_W,[+1,00,00],[00,+1,00],0,dim.ac_val,tp_rad,'φ'); # The point: output_point(op,dim,pos_W,dim.pt_rad,r_str,'θ','φ'); #---------------------------------------------------------------------- def output_system_SZ(op,dim,pos_W) : output_coord_axis(op,dim,[+1,00,00],[-5,+5],'A'); output_coord_axis(op,dim,[00,00,+1],[-5,-5],'Z'); if (op.rho) : r_str = 'ρ'; else : r_str = 'r'; # Relevant points: ooo_W = [00,00,00]; Pr_W = [00,00,dim.rc_val]; Pz_W = [00,00,dim.zc_val]; Pre_W = vec_from_ang_rad([00,00,+1],[+1,00,00],dim.ec_val,dim.rc_val); Qr_W = [dim.rc_val*sin(dim.ec_val),00,00]; Qra_W = [pos_W[0],pos_W[1],00]; Sr_W = [dim.rc_val,00,00]; Sra_W = vec_from_ang_rad([+1,00,00],[00,+1,00],dim.ac_val,dim.rc_val); # Radial traces: output_straight_trace(op,dim,ooo_W,Pr_W,True,r_str); output_straight_trace(op,dim,ooo_W,Pre_W,True,''); output_straight_trace(op,dim,ooo_W,pos_W,True,''); output_straight_trace(op,dim,Pz_W,Pre_W,True,''); output_straight_trace(op,dim,Pz_W,pos_W,True,''); # output_straight_trace(op,dim,ooo_W,Qr_W,True,''); # output_straight_trace(op,dim,ooo_W,Qra_W,True,''); output_straight_trace(op,dim,ooo_W,Sr_W,True,''); output_straight_trace(op,dim,ooo_W,Sra_W,True,''); # Vertical traces: # output_straight_trace(op,dim,Qr_W,Pre_W,True,''); # output_straight_trace(op,dim,Qra_W,pos_W,True,''); # Elevation traces: xy_W = vec_from_ang_rad([+1,00,00],[00,+1,00],dim.ac_val,1.0); output_circular_trace(op,dim,ooo_W,[00,00,+1],[+1,00,00],0,dim.ec_val,dim.rc_val,'θ'); output_circular_trace(op,dim,ooo_W,[00,00,+1],[+1,00,00],dim.ec_val,math.pi/2,dim.rc_val,''); output_circular_trace(op,dim,ooo_W,[00,00,+1],xy_W,0,dim.ec_val,dim.rc_val,''); output_circular_trace(op,dim,ooo_W,[00,00,+1],xy_W,dim.ec_val,math.pi/2,dim.rc_val,''); # Azimuth traces: tp_rad = dim.rc_val*sin(dim.ec_val); output_circular_trace(op,dim,Pz_W,[+1,00,00],[00,+1,00],0,dim.ac_val,tp_rad,'φ'); # output_circular_trace(op,dim,ooo_W,[+1,00,00],[00,+1,00],0,dim.ac_val,tp_rad,''); output_circular_trace(op,dim,ooo_W,[+1,00,00],[00,+1,00],0,dim.ac_val,dim.rc_val,''); # The point: output_point(op,dim,pos_W,dim.pt_rad,r_str,'θ','φ'); #---------------------------------------------------------------------- def output_coord_axis(op,dim,ax_dir,lab_dp,lab_str) : "Draws a coordinate axis." \ "\n" \ " The axis goes from the origin to {dim.ax_len*ax_dir}. The arrow tip has" \ " size {dim.hd_len}. The label is {lab_str} offset by {lab_dp} from the tip." # World points: ooo_W = [00,00,00]; # Origin. tip_W = rn.add(ooo_W,rn.scale(dim.ax_len,ax_dir)); # Tip of axis. # Image points: ooo = img_point(ooo_W,dim.map); # Origin. tip = img_point(tip_W,dim.map); # Axis tip. # Axis line: sys.stdout.write( \ ' <!-- The '+ lab_str +' axis: -->\n' \ ' <path d="M '+ pfm(ooo) +' L '+ pfm(tip) +'" stroke="rgb(0,0,0)"/>\n' \ ); # Tics: # Axis line: sys.stdout.write( \ ' <g stroke="rgb(0,0,0)" fill="rgb(0,0,0)">\n' \ ); ct = dim.ax_unit; while (ct < dim.ax_len - dim.hd_len) : tic_W = rn.add(ooo_W,rn.scale(ct,ax_dir)); # Axis tic (World). tic = img_point(tic_W,dim.map); # Axis tic (Image). sys.stdout.write( \ ' <circle '+ xyfm(tic,'c') +' r="1px"/>\n' \ ); ct += dim.ax_unit; sys.stdout.write( \ ' </g>\n' \ ); # Auxiliary directions perpendicular to the axis: bx_dir = [ax_dir[2], ax_dir[0], ax_dir[1]]; # Sideways direction B for arrow tip. cx_dir = [ax_dir[1], ax_dir[2], ax_dir[0]]; # Sideways direction C for arrow tip. output_arrow_head(op,dim,tip_W,ax_dir,bx_dir,cx_dir); lab_pos = rn.add(lab_dp, tip); output_label(op,dim,tip,lab_dp,make_italic_style(lab_str)); sys.stdout.write('\n'); #---------------------------------------------------------------------- def output_arrow_head(op,dim,tip_W,ax_dir,bx_dir,cx_dir) : "Draws an arrow head with tip at point {tip}, direction {ax_dir}, length {hd_len}." \ "\n" \ " The directions {bx,dir,cx_dir} must be orthogonal to {ax_dir} and to each other." # World points: eip_W = rn.add(tip_W,rn.scale(-dim.hd_len,ax_dir)); # Base of arrow head. tbm_W = rn.add(eip_W,rn.scale(-0.25*dim.hd_len,bx_dir)); # Corner B+ of arrow head. tbp_W = rn.add(eip_W,rn.scale(+0.25*dim.hd_len,bx_dir)); # Corner B- of arrow head. tcm_W = rn.add(eip_W,rn.scale(-0.25*dim.hd_len,cx_dir)); # Corner C+ of arrow head. tcp_W = rn.add(eip_W,rn.scale(+0.25*dim.hd_len,cx_dir)); # Corner C- of arrow head. # Image points: tip = img_point(tip_W,dim.map); # Arrow tip. tbm = img_point(tbm_W,dim.map); # Corner B+ of arrow head. tbp = img_point(tbp_W,dim.map); # Corner B- of arrow head. tcm = img_point(tcm_W,dim.map); # Corner C+ of arrow head. tcp = img_point(tcp_W,dim.map); # Corner C- of arrow head. sys.stdout.write( \ ' <path d="M '+ pfm(tcm) +' L '+ pfm(tbm) +' L '+ pfm(tcp) +' L '+ pfm(tbp) +'" stroke="rgb(0,0,0)" fill="rgb(0,0,0)" />\n' \ ' <path d="M '+ pfm(tip) +' L '+ pfm(tbm) +' L '+ pfm(tbp) +'" stroke="rgb(0,0,0)" fill="rgb(0,0,0)" />\n' \ ' <path d="M '+ pfm(tip) +' L '+ pfm(tcm) +' L '+ pfm(tcp) +'" stroke="rgb(0,0,0)" fill="rgb(0,0,0)" />\n' \ ); #---------------------------------------------------------------------- def output_label(op,dim,pos,dp,str) : "Ouputs a label {str} at image position {pos} displaced by {dp}." if (str != '') : dpx = dp[0]; dpy = dp[1]; # Compute anchor {anch}: if (dpx < 0) : anch = 'end' elif (dpx > 0) : anch = 'start' else : anch = 'middle'; if (dpy > 0) : # Adjust {dpy} for font height: dpy += 0.6*dim.font_wy; elif (dpy == 0) : dpy += 0.3*dim.font_wy; pos = rn.add(pos, [dpx,dpy]); sys.stdout.write( \ ' <text '+ xyfm(pos,'') +' stroke="rgb(255,255,255)"' \ ' stroke-width="3px" fill="none" text-anchor="'+ anch +'">'+ str +'</text>\n' \ ' <text '+ xyfm(pos,'') +' stroke="none"' \ ' fill="rgb(0,0,0)" text-anchor="'+ anch +'">'+ str +'</text>\n' \ ); #---------------------------------------------------------------------- def output_circular_trace(op,dim,ctr_W,u_W,v_W,ang_ini,ang_fin,rad,lab_str) : "Draws an angular coordinate trace.\n" \ "\n" \ " The arc goes from {d_ini} to {d_fin} given the center {ctr_W}, the plane's ortho" \ " axes {u_W,v_W}, the angle interval {ang_ini,ang_fin}, and the radius {rad}." \ "\n" \ " If the label is not empty, uses solid colored line, else a dashed " # Parameters: dashed = (lab_str == ''); ang = ang_ini; # Current angle. step = math.pi/100; # Initial guess for step. if (dashed) : down = False; # Is the pen down? dlen = 2; # Length of current dash/gap (px). color = '0,0,0'; else : down = True; # Is the pen down? dlen = 4; # Length of current dash/gap (px). color = dim.ct_color; # Compute the dash start and stop points {dini[k],dfin[k]}: dini = []; dfin = []; nd = 0; # Initial point: pos_W = rn.add(ctr_W, vec_from_ang_rad(u_W,v_W,ang_ini,rad)); # Current World point. pos = img_point(pos_W,dim.map); # Current Image point. ctr = img_point(ctr_W,dim.map); # Image arc center. # Loop on steps: while (ang < ang_fin) : pos_old = pos; # Save previous position. ang_old = ang; # Save previous angle. # Find {ang > ang_old} such that the step takes us about {dlen} avawy from {pre} while (True) : ang = ang_old + step; pos_W = rn.add(ctr_W, vec_from_ang_rad(u_W,v_W,ang,rad)); # Current World point. pos = img_point(pos_W,dim.map); # Current Image point. dst = rn.dist(pos_old,pos); rel = dst/dlen; # sys.stderr.write("ang = %6.2f dst = %6.2f dlen = %6.2f rel = %8.4f\n" % (ang,dst,dlen,rel)); if ((rel >= 0.9) and (rel <= 1.1)) : break; elif (rel < 0.5) : step *= 2; elif (rel > 2.0) : step /= 2; else : step /= rel; # Trim the angle to {ang_fin}: if (ang > ang_fin) : ang = ang_fin; pos_W = rn.add(ctr_W, vec_from_ang_rad(u_W,v_W,ang,rad)); # Current World point. pos = img_point(pos_W,dim.map); # Current Image point. # Process a gap/dash from {pre} to {pos}: if (down) : # Dash: dini[nd:nd] = [pos_old]; dfin[nd:nd] = [pos]; nd += 1; if (dashed) : down = False; dlen = 3; else : # Gap: down = True; dlen = 4; nd = len(dini); if (not dashed) : # Paint the sectors: sys.stdout.write( \ ' <g stroke="none" fill="rgb(100,0,200)" fill-opacity="0.25">\n' \ ); for k in range(nd) : sys.stdout.write( \ ' <path d="M '+ pfm(dini[k]) +' L '+ pfm(dfin[k]) +' L '+ pfm(ctr) +'"/>\n' \ ); sys.stdout.write( \ ' </g>\n' \ ); # Draw the arc: sys.stdout.write( \ ' <g stroke="rgb(' + color + ')" fill="none">\n' \ ); for k in range(nd) : sys.stdout.write( \ ' <path d="M '+ pfm(dini[k]) +' L '+ pfm(dfin[k]) +'" />\n' \ ); sys.stdout.write( \ ' </g>\n' \ ); # Compute the label position and displacement: ang_mid = (ang_ini + ang_fin)/2; mid_W = rn.add(ctr_W, vec_from_ang_rad(u_W,v_W,ang_mid,rad)); # Arc midpoint. mid = img_point(mid_W,dim.map); # Midpoint of arc. dmd,emd = rn.dir(rn.sub(mid,ctr)); # Direction from center to midpoint of arc. dpx = 5*dmd[0]; dpy = 5*dmd[1]; output_label(op,dim,mid,[dpx,dpy],make_italic_style(lab_str)); #---------------------------------------------------------------------- def output_straight_trace(op,dim,ini_W,fin_W,draw,lab_str) : "Draws and/or labels a dashed straight line.\n" \ "\n" \ " The line goes from {ini_W} to {fin_W} and is labeled {lab-str}. If {draw}" \ " is FALSE, does not plot it, just writes the label." \ "\n" \ " " # Parameters: dashed = (lab_str == ''); if (dashed) : style = 'stroke="rgb(0,0,0)" stroke-dasharray="3,4" stroke-dashoffset="9"'; else : style = 'stroke="rgb(' + dim.ct_color + ')"'; ini = img_point(ini_W,dim.map); # Start point. fin = img_point(fin_W,dim.map); # End point. if (draw) : sys.stdout.write( \ ' <path d="M '+ pfm(ini) +' L '+ pfm(fin) +'" ' + style + '/> -->\n' \ ); # Compute the label position and displacement: mid = rn.scale(0.5,rn.add(ini,fin)); # Midpoint of trace. dtr,ntr = rn.dir(rn.sub(ini,fin)); # Direction and length of trace. dpx = +5*dtr[1]; dpy = -5*dtr[0]; ooo = img_point([00,00,00],dim.map); # Origin. if (abs(dpx) > abs(dpy)) : if ((mid[0] < ooo[0]) != (dpx < 0)) : dpx = -dpx; dpy = -dpy; else : if ((mid[1] < ooo[1]) != (dpy < 0)) : dpx = -dpx; dpy = -dpy; output_label(op,dim,mid,[dpx,dpy],make_italic_style(lab_str)); #---------------------------------------------------------------------- def output_point(op,dim,pos_W,rad_W,x_str,y_str,z_str) : "Draws a dot with the given World position and radius, and its coordinate labels." # Image point and radius: pos = img_point(pos_W,dim.map); # Image point. rad = img_scale(pos_W,dim.map)*rad_W; # Image radius. if (rad > 0) : sys.stdout.write( \ ' <circle '+ xyfm(pos,'c') +' r="'+ dts(rad) +'"' \ ' stroke="rgb(0,0,0)" fill="rgb(0,0,0)"/> -->\n' \ ); # Compute the label position and displacement: ooo = img_point([00,00,00],dim.map); # Origin. if (pos[0] > ooo[0]) : dpx = +20 else : dpx = -20; dmd,emd = rn.dir(rn.sub(pos,ooo)); # Direction from center to midpoint of arc. dpy = 0; lab_fmt = make_roman_style('(') \ + make_italic_style(x_str) \ + make_roman_style(',') \ + make_italic_style(y_str) \ + make_roman_style(',') \ + make_italic_style(z_str) \ + make_roman_style(')'); output_label(op,dim,pos,[dpx,dpy],lab_fmt); #---------------------------------------------------------------------- def output_figure_postamble(op,dim) : "Writes the SVG postamble to {stdout}." sys.stdout.write( \ ' </g>\n' \ '</svg>\n' \ ) def vec_from_ang_rad(u,v,ang,rad) : "Converts polar coords {ang,rad} to Cartesia, given two axis directions {u,v} in {R^3}." c = rad*cos(ang); s = rad*sin(ang); return rn.add(rn.scale(c,u),rn.scale(s,v)) #---------------------------------------------------------------------- def dts(x) : "Converts the decimal number {x} to string." return ("%r" % x) #---------------------------------------------------------------------- def img_point(p,map): "Computes a two-dimensional Cartesian Image point {p} from a World point." if (len(p) == 3) : p = copy.copy(p); p[0:0] = [1]; assert (len(p) == 4), "invalid point"; q = rmxn.map_row(p,map); return [q[1]/q[0], q[2]/q[0]]; def img_scale(p,map): "Computes the World to Image magnificatin factor at the World point {p}." if (len(p) == 3) : p = copy.copy(p); p[0:0] = [1]; assert (len(p) == 4), "invalid point"; q = rmxn.map_row(p,map); s = sqrt(map[1][1]*map[1][1] + map[2][1]*map[2][1] + map[3][1]*map[3][1]) return s*p[0]/q[0]; def make_italic_style(str) : "Packages a string with italic style markup." if (str != '') : return '<tspan font-style="italic">' + str + '</tspan>' else : return '' def make_roman_style(str) : "Packages a string with normal style (non-italic) markup." if (str != '') : return '<tspan>' + str + '</tspan>' else : return '' def pfm(p) : "Formats a two-dimensional Image point {p} as two numbers separated by a comma." return "%+06.1f,%+06.1f" % (p[0],p[1]); def xyfm(p,tag) : "Formats a two-dimensional Image point {p} as '{tag}x=\"{p[0]}\" {tag}y=\"{p[1]}\"'." return "%sx=\"%+06.1f\" %sy=\"%+06.1f\"" % (tag,p[0],tag,p[1]); #---------------------------------------------------------------------- # Main program: op = parse_args(pp); output_figure(op);