/DFT_geom module | toolbox |============================== DFT_geom ================================= | | NOTE: Grid nodes are located in the center of their axial length interval. | Variables are thought to change linearly between grid nodes (or | values at grid nodes constitute the average of their interval). | | Boundary coordinates of grid subdomains must be integer multiples | of the smallest interval used in the grid. | | front end procedures: | | - build_pab_geom (pore/atrium/bulk) | 200 dict dup begin | DFT_geom |===================== Build pore/atrium/bulk geometry ================== | | - initializes the geometry arrays, X, Radius, Area, Height, and fills in | their values, thus defining a system geometry | | - invoke in a 'stopped' context (inconsistent geometry specifications) | | - expects the following the values of the following items in the channel | dictionary: | | Xslice - axial interval (in pore proper) | Poreradius - the radius of the cylindrical pore proper | Lengths - array of geometrical lengths: left bulk, | left atrium, pore proper, right atrium, right bulk | comp_L_angle - a procedure that computes the left atrial angles | comp_R_angle - a procedure that computes the right atrial angles | comp_P_radius - a procedure that computes the radii in the pore | | 'Angle' procedures map angles from 0 to Pi/4 to the X values 0 to 1. | The 'radius' procedure maps pore radii (in m) on the X values 0 to 1. | | - the system includes two hemispherical bulk zones interfacing with a | (cylindrical or arbitrarily shaped) pore proper via two funnel-like | atrial domains | - each zone is described as a spherical cone of appropriate aperture | angle | - the grid abscissa corresponds to the axis intersections of the | spherical cones surfaces | - the grid comprises a uniformly spaced central pore section | and two external sections whose spacing is increased toward the | periphery in proportion to system area | | Creates the following information in the channel dictionary: | | SystemLength - overall axial length of system | cLengths - of cone segments on pore axis | X - normalized (to SystemLength) axial coordinates of grid | DX - slice widths of grid | H - normalized increments of axial coordinate | Radius - pore radius (defined over atrial sections and pore | proper) | XP - abscissae to which the Radius applies | Area - equipotential surface area of the system that intersects | the axis at the corresponding grid node (defined over | full system); equipotential surfaces are designed to | intersect the pore wall at a right angle | Indices - indices of the first grid nodes in the sections described | in cLengths | Dimensions - number of nodes in each section described in cLengths | poreproper - carving operator: carves out poreproper from grid array | leftabath - carving operator: left atrium and bath | rightabath - carving operator: right atrium and bath | /build_pab_geom { |-- build temporaries and compute geometry arrays save /temp_dft_var name /poreradius Lengths 2 get inslices /d array def poreradius channel comp_P_radius _channel /Poreradius poreradius 0 get def /phi Lengths 1 get inslices /d array def phi makecenters comp_L_angle pop Lengths 0 get inslices Lengths 1 get inslices add comp_tabs comp_surface /LXP name /LA name /LR name /LX name /LDX name /Poreradius poreradius last get def /phi Lengths 3 get inslices /d array def phi makecenters comp_R_angle pop Lengths 3 get inslices Lengths 4 get inslices add comp_tabs comp_surface /RXP name /RA name /RR name /RX name /RDX name /aux Lengths 2 get inslices /d array def | for pore proper temp_dft_var capsave define_final_grid channel DX _channel 0 LDX flip fax Xslice aux copy fax RDX fax pop pop DX WA copy integrateRS pop /cLengths [ WA Indices 1 get 1 sub get WA Indices 2 get 1 sub get WA Indices 1 get 1 sub get sub WA Indices 3 get 1 sub get WA Indices 2 get 1 sub get sub WA Indices 4 get 1 sub get WA Indices 3 get 1 sub get sub WA last get WA Indices 4 get 1 sub get sub ] channel def _channel channel X _channel 0 LX flip cLengths 0 get cLengths 1 get add /xoffs name neg xoffs add fax aux makecenters cLengths 2 get mul xoffs add fax /xoffs xoffs cLengths 2 get add def RX xoffs add fax pop pop channel XP _channel 0 LXP flip Lengths 0 get Lengths 1 get add /xoffs name neg xoffs add fax aux makecenters Lengths 2 get mul xoffs add fax /xoffs xoffs Lengths 2 get add def RXP xoffs add fax pop pop channel Area _channel 0 LA flip fax poreradius aux copy dup mul Pi mul fax RA fax pop pop channel Radius _channel 0 LR flip fax poreradius fax RR fax pop pop | temp_dft_var restore channel define_XH define_zones _channel } bind def |-------------------- internal procedures of geometry builders ----------------- |----- place ramp of interval centers (normalized) into the given array | arr | arr /makecenters { /dest name /ndest dest length def /step 1.0 ndest div def /first step 0.5 mul def dest 0 ndest first step ramp pop } bind def |----- place ramp of interval right boundaries (normalized) into the given array | arr | arr /makerightbs { /dest name /ndest dest length def /step 1.0 ndest div def /first step def dest 0 ndest first step ramp pop } bind def |----- define X and H profiles /define_XH { 0.0 Lengths add /SystemLength name XP SystemLength div pop X SystemLength div H copy 1 H length 1 sub getinterval X 0 X length 1 sub getinterval sub pop } bind def |----- define final grid and some profile buffer arrays /define_final_grid { channel /Ngrid Lengths 2 get inslices LA length add RA length add def /NLB LA length LR length sub def /NRB RA length RR length sub def /Indices [ 0 NLB LX length dup Lengths 2 get inslices add dup RX length NRB sub add ] def /Dimensions [ NLB LX length NLB sub Lengths 2 get inslices RX length NRB sub NRB ] def [ /XP /Radius ] { Ngrid NLB sub NRB sub /d array def } forall [ /Area /Height ] { Ngrid /d array def } forall [ /DX /X /H /WA /WB /WC /WD ] { Ngrid /d array def } forall _channel } bind def |----- define zone operators /define_zones { /poreproper [ Indices 2 get Dimensions 2 get /getinterval mkact ] mkact bind def /leftabath [ Indices 0 get Dimensions 0 get Dimensions 1 get add /getinterval mkact ] mkact bind def /rightabath [ Indices 3 get Dimensions 3 get Dimensions 4 get add /getinterval mkact ] mkact bind def } bind def |----- reverse order of elements in an array | arr | rev_arr /flip { /flip_arr name /flip_hi flip_arr length 1 sub def 0 1 flip_arr length 2 div 1 sub { /flip_lo name flip_arr flip_hi get flip_arr flip_lo get flip_arr flip_hi put flip_arr flip_lo put /flip_hi flip_hi 1 sub def } for flip_arr } bind def |----- compute an integer dimension from a real range /inslices { Xslice div 1e-8 add /l ctype } def |----- compute the radius and surface in a conical domain | requires the phi (angles) array, creates its own output arrays | | Some compromises are necessary to avoid conflicts in the design of | grid intervals: (1) the projection of xp to x in the atrial zone can | in principle exceed the increments of x that are demanded by the | increase in area in the adjacent bulk; we force bulk intervals to be | no smaller than the last atrial interval; (2) if interval boundaries | are to exactly match the atrial/bulk boundary, a small grid interval can | result (so grid intervals are not monotonic); again, we force the atrial | boundary to move out such that the grid interval does not contract. Thus | the overriding rule is that grid intervals increase monotonically from | the pore proper outwards. A result is that the atrial boundary as defined | by Indices, Dimensions can be less than the limit set in Lengths. /comp_surface { /dx nab /d array def /a0 Pi Poreradius dup mul mul def /x_1 0.0 def /xp_1 0.0 def /area_1 a0 def /radius_1 Poreradius def /kx 0 def /kp 0 def /kv 0 def /slice Xslice def |-- up to atrial boundary as defined by spherical cone surface { { x kx get x_1 sub dx kv put x kx get dup x_1 add 0.5 mul x kv put /x_1 name xp kx get dup xp_1 add 0.5 mul xp kv put /xp_1 name radius kx get dup radius_1 add 0.5 mul radius kv put /radius_1 name area kx get dup area_1 add 0.5 mul area kv put /area_1 name /oldslice dx kv get def /slice area kv get a0 div Xslice mul dup oldslice lt { pop oldslice } if def /kv kv 1 add def /kx kx 1 add def { kx na 1 sub ge { stop } if x kx get x_1 sub slice ge { exit } if /kx kx 1 add def } loop } loop } stopped pop /kp kv 1 sub def |-- up to end of bath { { { x kx 1 add get x_1 sub slice ge { exit } if /kx kx 1 add def kx nab 1 sub ge { stop } if } loop x kx get x_1 sub dx kv put x kx get dup x_1 add 0.5 mul x kv put /x_1 name area kx get dup area_1 add 0.5 mul area kv put /area_1 name /oldslice dx kv get def /slice area kv get a0 div Xslice mul dup oldslice lt { pop oldslice } if def /kv kv 1 add def /kx kx 1 add def } loop } stopped pop x kx get x_1 sub dx kv put x kx get x_1 add 0.5 mul x kv put area kx get area_1 add 0.5 mul area kv put /kv kv 1 add def dx 0 kv getinterval x 0 kv getinterval radius 0 kp getinterval area 0 kv getinterval xp 0 kp getinterval } bind def |-- compute radius, area, hs and rp as function of xp /comp_tabs { /nab name /na phi length def /xp na /d array def /x nab /d array def /hs na /d array def /area nab /d array def /radius na /d array def /rp na /d array def x makerightbs Xslice x length mul mul 0 na getinterval xp copy pop phi radius copy tan integrateRS Xslice mul Poreradius add pop phi rp copy sin -1 pwr radius mul pop phi hs copy cos neg 1.0 add rp mul pop rp area copy hs mul Pi 2 mul mul pop /xpap x na 1 sub get def x 0 na getinterval hs add pop /xpa x na 1 sub get def false na 1 nab 1 sub { /k name x k get xpa gt { pop true exit } if } for { /n nab k sub def x k n getinterval x na n getinterval copy pop x na n getinterval area na n getinterval copy xpap sub dup mul 2.0 Pi mul mul pop /nab na n add def } { (\n** Bulk boundary inside atrial cone! **) toconsole stop } ifelse } bind def |------------------------ some front-end tools --------------------------- |----------------------------- smear a number over the volume of a segment | The result is a volume density | | amount destination first dimension | -- /smear_volume { [ 3 1 roll /getinterval mkact ] mkact /seg name /smear_dest name /amt name Area seg smear_dest seg copy DX seg mul integrateRS last get /vol name amt vol div smear_dest seg copy pop } bind def end _module