4 x.- 22):: ..F.t:1. i» P! .., 5 a .1 J h 3.3:); 3.3.... r {23.1. n, . . x‘lf... IVERSITY LIBRARIES llljllllflllilflm u m: m u 897 7666 This is to certify that the dissertation entitled Probabilistic Models to Represent Loads Due To Human Activities presented by Bassem K. Khafagi has been accepted towards fulfillment of the requirements for Doctor of Philosophy degree in Civil Engineering Major professor Mew MSU i: an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY Michigan State Unlverslty PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE MSU Is An Affirmative Action/Equal Opportunity Institution 7 , , ., c:\cilt\dBMLB.Dm3-DJ PROBABILISTIC MODELS TO REPRESENT LOADS DUE TO HUMAN ACTIVITIES By Bassem K. Khafagi A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Civil and Environmental Engineering April 1993 ABSTRACT PROBABILISTIC MODELS TO REPRESENT LOADS DUE TO HUMAN ACTIVITIES BY Bassem K. Khafagi Human loads are random and produce dynamic force components which cannot be predicted with certainty. Current codes account for human loads by assuming a static load equivalence that included the maximum dynamic effect of human movements. The primary objective of this research is to generate probabilistic functions which accurately represent the dynamic force caused by several human actions. Measured data were obtained from three sets of experiments, two designed to measure forces due to the action of an individual on a force platform; and the third designed to measure the forces due to the action of a group on a floor system. The study included six different in-place human movements grouped into two categories; transient motions (single jump, standing up suddenly, and dropping into a seat), and periodic motions (jumping, jouncing, and swaying side to side). Transient actions were modeled by passing straight or curved lines through control points that define the force-time history. Best-fit distributions for the control points were determined. Model error was developed as a Gaussian random process and was incorporated into the model to better represent the modeled actions. Periodic loadings were analyzed as a series of cycles where each cycle is a full period of the motion considered. The period and peak force ratio of each cycle were modeled as first order autoregressive processes. Model verification was applied to the different motions considered in the study. Stochastic models were compared with experimental data and feedback from the process was used to refine the models. A computer program; Human Load Simulation, HLS, was developed to simplify generating load distributions for any combination of human loading conditions. Output of the program includes the statistical parameters of the control points associated with the time histories, amplitudes of peaks and their arrival time, spectral energy distributions for periodic motions, and force—time history for each simulation. A modified version of HLS was developed to simulate group loading conditions. Several potential applications of HLS were presented. Code values were examined and found to be acceptable except for the transverse force component where a value of 45 lb/ft was recommended to replace the current 10 lb/ft value in the code. An example showing the advantage of using probabilistic loads was presented. Other potential applications of the research outcome were discussed. It is anticipated that research findings will have an impact on standards for structural design and serviceability control. Copyright" by BASSEM KAMAL KHAFAGI April 1993 To Ekram, Yasmeen, Louai, and ........ iv ACKNOWLEDGE/IENTS In the course of this project, I have received assistance in a variety of forms from many pe0ple. I am indebted to them, and I wish to acknowledge this debt. My heartfelt thanks to Professor W. Saul and Professor R. Harichandran for their cooperation, enlightening advice and devoting part of their time to guide and help me to overcome the difficulties I have faced during my research program; Professor T. Wolf, Professor F. Hatfiled, and Professor G. Ludden, other members of my Graduate Committee, for their recommendations and review of the dissertation; Special thanks go to Ahmad Al- Ghamedi and Abdulrahman Al-Hadlag for their friendship and incessant help. I am grateful to my enlightened parents who supported and encouraged me throughout my education. Finally, I believe that God is behind any success I have reached. So, my thanks to Him. ..... TABLE OF CONTENTS LIST OF TABLES ..................................... viii LIST OF FIGURES ....................................... x CHAPTER 1.0 INTRODUCTION ............................. 1 1.1 General ...................................... 1 1.2 Problem Statement ................................ 2 1.3 Objectives of the Study ............................. 3 1.4 Research Description ............................... 5 CHAPTER 2.0 CURRENT STATE OF KNOWLEDGE ............... 8 2.1 Introduction .................................... 8 2.2 Literature Review ................................ 8 2.3 Current Project ................................. 14 2.4 Current Practices ................................ 16 CHAPTER 3.0 EXPERIMENTAL MEASUREMENTS .............. 18 3.1 Instrumentation ................................. 18 3.2 Data Collection ................................. 21 3.3 Data Processing ................................ 24 3.4 Statistical Output Summary .......................... 27 CHAPTER 4.0 MEASUREMENT ANALYSIS .................... 35 4.1 General ..................................... 35 4.2 Methodology .................................. 35 4.3 Transient Loading ................................ 39 4.3.1 Single Jump .............................. 40 4.3.2 Sitting Down ............................. 50 4.3.3 Standing Up ............................. 53 4.4 Periodic Loading ................................ 55 4.4.1 Periodic Jumping .......................... 57 4.4.2 Periodic Jouncing .......................... 63 4.4.3 Side Swaying ............................. 64 vi CHAPTER 5.0 HUMAN LOAD MODELING .................... 67 5.1 Introduction ................................... 67 68 ................................. 5 .2 Basic Concepts 5.2.1 Random Processes .......................... 68 5.2.2 Random Number Generation ................... 71 5 .2.3 Generation of Correlated Random Variables .......... 73 5.2.4 Monte Carlo Simulation ...................... 74 5.3 Transient Load Modeling ........................... 75 5.3.1 Mathematical Approach ...................... 76 5.3.2 Load Simulation ........................... 85 5 .4 Periodic Load Modeling ............................ 86 5.4.1 Stochastic Analysis ......................... 88 5.4.2 Load Modeling ............................ 92 CHAPTER 6.0 MODEL APPLICATIONS ...................... 95 6.1 Introduction ................................... 95 6.2 Group Simulation ................................ 95 6.3 Code Verification ................................ 99 6.4 Reliability-Based Design ........................... 102 6.5 Other Potential Applications ........................ 105 CHAPTER 7.0 CONCLUSIONS ............................ 107 7.1 Summary .................................... 107 7.2 Conclusions .................................. 109 7.3 Future Work .................................. 110 APPENDIX A DATA PROCESSING PROGRAMS ................ 112 APPENDIX B LOAD SIMULATION PROGRAMS ......... . ...... 159 LIST OF REFERENCES ................................. 186 vii Table 2.1 2.2 2.3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 Lainipioi-u LIST OF TABLES Periodic actions .................................. Representative activity and their applicability to actual structures ..... Bounds for Natural Frequencies of Structures subject to Human Loads . . Transient Data-base ................................. Periodic Actions Data-base ............................ Typical Output Files from "DATA.FOR" ................... Output summary of Single Jump (One~Individual) .............. Maximum Force Ratios (Single Jump) ..................... Peak Statistics for Floor System Data ...................... 2 Hz.—Periodic Jumping Peak Analysis ..................... Human Movements Considered in The Study ................. Human Movements Considered in The Study ................. Control Points for Single Jump .......................... Single Jump Control Points Statistics ...................... Correlation Matrix (Single Jump) ........................ Best-Fit Distributions (Single Jump) ....................... Sitting Down Control Points Statistics ...................... Correlation Matrix (Sitting Down) ........................ Best-Fit Distribution (Sitting Down) ....................... Standing Up Control Points Statistics ...................... Best-Fit Distribution (Standing Up) ....................... Correlation Matrix (Standing up) ........................ 2-Hz. Periodic Jumping .............................. Peak Analysis (Periodic Jouncing) ........................ Control Points Statistics (Periodic Jumping) .................. Periodic Jouncing Statistics (One Individual, 2 Hz.) ............. Max. Peaks (Side Swaying) ............................ Typical Side Swaying Data ............................ Measured Versus Simulated Control Points for Single Jump ........ Measured Versus Simulated Control Points For Sitting Down ....... Measured Versus Simulated Control Points For Standing up ........ Periods and Peaks for a Typical Periodic Jumping .............. Constants a1, and aT for the Periodic Motions ................. viii 14 17 22 33 35 37 44 46 47 5 1 52 53 54 54 55 57 61 63 63 66 67 67 87 88 89 94 95 Table Page 5 .6 Measured Versus Simulated Data (Periodic Jumping) ............. 95 6.1 Measured Versus Simulated Maximum Force Ratios ............. 98 6.2 Measured Versus Modified Simulated Maximum Force Ratios ....... 100 6.3 Dynamic Force Ratios for One Individual ................... 104 Figure 1.1 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 LIST OF FIGURES Page Human Loads As Static and Dynamic Forces .................. 2 View of the Force Platform, Saul et. a1. ( 1990) ................ 19 The Floor System Lab, Saul et.al. (1990) .................... 20 Human Load Experimental Data-base ...................... 21 DATA.FOR Flow Chart .............................. 24 Force Type History (DATA.FOR Screen Output) ............... 27 Spectral Analysis (DATA.FOR Screen Output) ................ 27 A Typical Force-Time History (Floor.For Output) .............. 31 Mean, Min, Max Peaks For 2 Hz. Periodic Jumping ............. 33 Log-Linear Regression Model ........................... 34 A Single Jump .................................... 38 Control Points for A Single Jump ......................... 39 Control Points for a Periodic Motion ....................... 39 A Typical Single Jump ............................... 41 Best-Fit Distribution for Response Lag ..................... 50 A Typical Sitting Down Motion .......................... 53 A Typical Standing up ............................... 55 A Typical Periodic Jumping ........................... 59 Control Points for a Typical Periodic Jumping ................ 6O CONTROLFOR Output (Periodic Jumping, 2 Hz. , Parallel Ratio) . . . . 62 CONTROLFOR Output (Periodic Jumping, 2 Hz., Vertical Ratio) . . . . 62 A Typical Periodic Jumping Response Period ................. 64 A Typical Periodic J ouncing ............................ 65 A Typical Swaying Side to Side .......................... 66 Typical Samples in the Random Process F(t) .................. 71 HLS Screen Output (Single Jump, One Individual) .............. 78 HLS Screen Output (Sitting Down, One Individual) .............. 78 Error Function (Single Jump, One Individual) ................. 82 Error Spectrum (Sitting Down, One Individual) ................ 82 Straight and Curved Lines of Single Jump Model ............... 83 Single Jump Error Function, b = 0.73 Hz. ................... 85 Single Jump Error Function, b = 0.36 Hz. ................... 85 Figure Page 6.1 A Typical GROUP.FOR Output, 10 Individuals ................ 99 6.2 A Typical GROUP.FOR Output, 40 Individuals ................ 99 xi 1.0 INTRODUCTION 1.1 General Human loads are random and therefore cannot be predicted with certainty. In addition to the forces due to static weights, the movement of persons produces dynamic force components. Human loads vary from event to event in an apparent pattern and therefore it was postulated they could be shown to obey the laws of probability and hence would be predictable. In addition to many other possible applications, these loads are the primary concern in the design and analysis of assembly structures. A challenge a structural engineer faces is how to quantify the uncertainty associated with loads due to , human movements and provide corresponding measures to assure safety and serviceability of the structure at a reasonable cost. This was the initial motivation for attempting to provide procedures for generating functions to represent forces due to human movements. When a person on a horizontal surface is absolutely motionless a vertical force is brought to bear on the support due to gravitational acceleration acting on the mass of the person's body, the weight of the person, as shown in Figure 1.1.a. The value of the static force could be predicted by a probability distribution model of the individual weights of a population sample of pe0p1e. When activity takes place, such as jumping, sitting down, or swaying, there are resultant vertical and horizontal forces which fluctuate ith time due to the changing acceleration of the person's mass, as illustrated in Figure 2 1.1.b. The vertical component of this force varies from zero to several times the person's weight. Horizontal components do not exist in the static situation. The time variation for any one of various human actions, for instance a single jump in place, provides a pattern characteristic of that action. Similar to the distribution model for No Force Force Time Time Force Component h the x-dlr Force Convonent h the x-dlr Face NoFOmQ Rm W Tune Tune ForceComponenthmey-dlr ForceComponenthuIey-dlr Time Time Force Component h the z-dk Force Component in the z-dlr Figure 1.1 Human Loads As Static and Dynamic Forces predicting the weight of an individual in a group of people, the parameters defining the characteristic pattern of a force due to a human action vary randomly within the bounds of distribution models. 1.2 Problem Statement For the design of structures to avoid catastrophic failure, the uncertainty of human oads was Traditionally accounted for by assuming a static load equivalence that included e maximum dynamic effect of human movements. Recent incidents of occupant 3 discomfort, and strong vibration of assembly structures, Bachmann (1992), have reinforced the need for quantification of dynamic human loadings. In recent years, researchers have been utilizing probabilistic methods to model loads and calculate the probabilities associated with combined load effects. A number of investigations have been conducted on human loads variabilities which have resulted in several important achievements, Allen (1990), Bachmann (1992), and Murray (1991). However, up until this study, no satisfactory method existed for determining the statistics of the dynamic load one person in a group would generate when performing various types of human movements. The focus of this study is the modeling and analysis of human loadings to include and quantify their random nature. The impact of the research on the evaluation of safety of assembly structures and determination of design loads and load combinations in building codes is summarized. 1.3 Objectives of the Study The primary objective of the research is to generate probabilistically-based functions which will accurately represent the dynamic force pattern and values caused by several types of human actions. The force due to each person in a group of n persons is of the form {F(t)} where F,(t) represents a single force and i = 1, 2, 3, . . . ,n. The parameters ceded to generate these probability-based-functions will be determined utilizing measured alues. Data from tests conducted by the author are supported with data from earlier tages of the research. The force functions produced will represent the dynamic load of 4 an individual in a crowd of people who are also active. With accurate correlations, the entire group action will then be modeled with a matrix of discrete forces. In addition to the interest in having these forces derived and available, there are a number of ancillary objectives which will have benefits and future possibilities: I The models will be especially important in cases where sustaining human loading is the primary purpose of the structure such as stadium or grandstand seating. In these cases the models would be valuable for design and analysis. I The models may be used to determine limits or maxima; for example, for the verification of specific building code values which are based on avoidance of catastrophic failure. The code values used in the design of grandstand seating areas, such as for stadium design, will be derived and compared with published code values for verification. I The models provide for the ability to conduct simulations to study any of a number of possible events. An example is the post-elastic response of a grandstand element subjected to a person jumping on it. I The models are valuable for reliability-based design methods where the design input needs to be probabilistically presented. I The load models may eventually have a major impact on standards for structural design such as those produced by ANSI, the American National Standards Institute, from which many building codes are derived. I The models will be necessary for the design of active controls to minimize vibration of structures subjected to human movements. Accurate definition of 5 time-varying forces where humans are the primary source of the forces, such as for the design of stadiums, are essential for providing active damping of the structure. I The models may be used to define parameters, such as acceleration, based on serviceability of the structure. Serviceability criteria should be a requirement for acceptable design. I It is expected that the process evolved in developing the load models will provide a methodology for extending the list of models to additional human actions such as for aerobics, climbing or descending stairs, and the " heel drop" test. I The models and the methodology should provide the insight needed to extend definitions of load functions to forces from actions due to persons moving such as dancing, walking or marching. 1.4 Research Description The load modeling of human movements was accomplished through four tasks: a) Experimental measurements: Force histories from tests performed on a small force platform and on a large floor system were used to statistically characterize dynamic loading due to human activities. b) Model identification: Digitally recorded responses due to specific motions were eviewed to obtain all available information on how the force histories are generated and hat factors affect the shape and magnitude of the generated force . m'fid” L L __-__r. 6 c) Load Modeling: The experimental data were analyzed and inferences were made regarding parameters needed to develop the load models. Modeling methodologies are presented along with detailed procedures for simulating the different human actions considered in the study. (1) Model verification : The proposed models were used to generate extensive load histories which were compared statistically to the measured test data. The intent of the process was to reveal inadequacies of the model and therefore, through iteration, improve the proposed model. This dissertation is organized in the following parts: Introduction (Chapter 1.0): It is a condensed outline of the topics herein, defined and more fully explained in the main body of the dissertation. It is intended as a broad mapping of areas covered by the research. Current State of Knowledge (Chapter 2.0): A detailed review of the current state of knowledge in the field of modeling human loads is presented. Earlier development of the project of which this research is part of is outlined. The chapter also includes an overview of current practice in the design process and how the present U. S. codes characterize human loads on assembly structures. Experimental Measurements (Chapter 3. 0): An overview of procedures of data ollection, the experimental setup, data manipulation, and data representations is resented. Computer programs were developed to automate the process of converting aw voltage data into time-force histories. Data available from earlier measurements ere utilized as well as data from new tests. 7 Measurement Analysis (Chapter 4. 0): For each human movement type chosen for the study, the forcing function model is hypothesized. The parameters needed to define each motion are evaluated and analyzed statistically. Probabilistic Load Modeling (Chapter 5.0): The chapter describes an iterative model- building methodology whereby the stochastic models introduced in Chapter 4.0 are related to actual time series data. The process of model verification was applied to the different motions considered in the study. A model to simulate group loading is provided. The developed forcing functions are compared with experimental data. Feedback from this process is used to refine the forcing functions. Model Applications (Chapter 6. 0): The impact of the research findings on the current state of knowledge is described. Several Applications showing the use of the Human Load Simulation program, HLS, are developed from the study. Recommendations on handling dynamic human loads in the design stage are presented. Conclusions (Chapter 7. 0): A summary of the research along with the major findings and conclusions is presented. The chapter also includes recommendations for designers and code agencies for incorporating research findings in future reliability-based-codes. A list of tasks for future research is outlined. . .~.---‘~.;v‘ mum- ' - 2.0 CURRENT STATE OF KNOWLEDGE 2.1 Introduction: In this chapter, a detailed review of the current state of knowledge in the field of modeling human loads is presented. First, an overview of the history of measuring force parameters of human loads is demonstrated. Then earlier deve10pment of the project of which this research is part of is outlined. Finally, the chapter outlines the current practice in both the design and analysis phases that relates to human load modeling. 2.2 Literature Review Studies of human live loads, from before the turn of the century, have been reported in depth by Saul and Tuan (1986). Tilden (1913) is the first person reported to have measured the parameters of forces due to human movements; his work formed the basis for the parameters used in present day building codes for grandstands and similar places of assembly. Much of the intervening work focused on finding the statically equivalent code values used for design today (ASA 1941). A series of incidents involving the collapse of portable grandstands in the early 1920's led to tests conducted under the supervision of the American Standards Association (now American National Standards Institute) to measure horizontal forces exerted by spectator movements. Other code 9 associations adopted the design values recommended by the American Standards Association Committee in 1941, and the values remain the same today. These tests focused on determining the upper bounds forces and neglected dynamic content. Bachmann (1984) examined the hypothesis that resonance behavior in a structure can be avoided if its natural frequency is higher than that of the forcing frequency and discovered that the assumption is not always correct. He found that lightly damped structures may be forced into significant resonance if their natural frequency is an integer multiple of the forcing frequency. A case study of a gymnasium where vibrations with large amplitudes due to resonance occurred at twice the forcing frequency is discussed in detail in his paper. He proved that the second and third harmonic of the forcing frequency may produce significant dynamic loading. Allen, et al (1975, 1982,1987) found that although structures in most cases were built to code, serviceability problems such as unacceptable vibrations due to occupants' movements occur. They related structural resonance to thirteen cases of serviceability problems caused by dancing, spectator activity at rock concerts, aerobics, exercising, and walking. Their studies were concerned with response of the structure and gave negligible treatment to the input function. Greimann and Klaiber (1978) generated a continuous forcing function to successfully model a resonance condition caused at Iowa State University stadium by spectators moving together on a grandstand. Tuan and Saul (1985) focused upon measuring the loads produced by an individual performing prescribed maneuvers on a small force platform. The loads produced by an individual performing prescribed maneuvers measured on a small force platform were 10 categorized and some statistical analyses were done on the data to produce tentative load functions. The research was continued by McDonald (1984) who designed and constructed a 4 by 8 foot force platform for measuring loads produced by up to five people. The platform was used to acquire data for actions due to one or more persons. Ebrahimpour (1989) , proposed group forcing functions to produce an equivalent vertical dynamic load to be compared to code presented values. A sinusoidal Fourier series was used to describe the total force developed by a symmetrically placed group of people jumping in unison at a fixed rate. VanKleek (1988), constructed and instrumented a 12 foot by 15 foot floor system to accommodate up to 40 people. The floor system was used to obtain experimental measurements and test the accuracy of the group load model pr0posed by Ebrahimpour. Statistical modeling was performed for periodic jumping conducted at frequencies of 2 and 3 Hz and also for a single jump. The results showed that for the specific conditions of the test, the equivalent maximum uniform floor load approached an asymptotic value comparable to the code value. Rainer, et a1 (1988), investigated the dynamic loading and response of footbridges using a force platform and asking individuals to walk or run along the span at a predetermined pulse rate. They concluded that for walking, running, and jumping, the excitation force can be represented as F(t): N F(t) = P [ 1 + E ansin(2mtfi + 91.)] .......... [2.1] n=l where: P = static weight of a person; a = Fourier amplitude or coefficient; 11 n = harmonic order of footstep rate; = footstep rate in steps per second; t = time in seconds; ¢ = relative phase angle; and N = total number of harmonics. Allen (1990) introduced some revisions to the Canadian Building Code design criteria Table 2.1 Periodic actions* 3.; .;.;. I Walking I Footbridges I Running , I Office buildings I Jumping I Gymnasia and sports halls I Dance and concert balls with no fixed I Dancing seating I Hand clapping with body bouncing while standing I Concert halls with fixed seating as well I Hand clapping only as galleries I Lateral body swaying while seated or while standing I High-diving platforms in swimming pools " Afler Bachmnn (1992). for floor structures to reduce building vibrations due to human activities. He studied a number of recent vibration problems in concrete building structures due to dancing and aerobics exercises and suggested the following design criteria for minimum natural frequency: f 2 f 1 + 1.3 . E22 ............... [2.2] O ao/g to where: = forcing frequency; a = natural frequency of the structures; l2 a/g = acceleration ratio; o ---= mass weight; and a», = Weight of the person. a = dynamic load factor He recommended increasing the factor (1 .3) to (2.0) in the case of aerobics activities. Considering the problem of floor vibrations, Murray (1991) provided a methodology for controlling annoying floor movement in residential, office, commercial, and gymnasium type environments. His criteria states that the motion of floor systems would not be objectionable to occupants of residential and office environments if the following inequality is met: D > 35 A0 f + 2.5 ................. [2.3] where D = damping in percent of critical, A, = maximum initial amplitude of the floor system due to heel-drop excitation (in.) defined as the force by a person standing on the toes and dropping suddenly from that position, and = first natural frequency of the floor system (Hz). For commercial environments (such as shopping centers), the criteria is based on an acceleration tolerance limit of 0.005 g due to walking excitation so that maximum deflection under 450 lbs. force applied anywhere on the floor system does not exceed 0.02 in. Murray (1991) confirmed the usefulness of Eq. [2.2] proposed by Alan (1990), to be used for gymnasium environments. Bachmann (1992) suggested the following function for the dynamic forces due to rhythmical body motions shown in Table 2.1: 13 Fp(t)=G+AGlsin(2rrj;t)+Astin(41tj;t- 4.8 Hz (C > 5%) and > 7.2 Hz (C < 5%)resp. Gymnasia and sports halls > 7.5 Hz > 8.0 Hz > 8.5 Hz > 9.0 Hz Dance and concert halls > 6.5 Hz > 7.0 Hz > 7.5 Hz > 8.0 Hz without fixed seating Concert halls and theaters Vertical: > 6.5 Hz with fixed seating (incl. Horizontal: > 2.5 Hz pop concerts) lateral sway bracing loads of 24 pounds per foot parallel to the seating and 10 pounds per foot perpendicular to seat and footboards for grandstands, reviewing stands, and bleachers. It also recommends using a uniform vertical live load of 100 pounds per square foot for assembly structures with moveable seating. A design criteria for assessing the acceptability of floor structures subjected to human periodic movements was introduced into the Commentary to the 1985 National Building Code of Canada. l7 Bachmann (1992) showed that rhythmical human body motions can induce considerable dynamic forces on assembly structures. Table 2.3 presents bounds for the natural frequency of structures subject to man-induced vibrations to reduce the probability of resonance when the forcing frequency approaches the natural frequency of a structure. In general, recent studies to include the dynamic effect of human loads represents an advance in representing the effect of human loads. However, work was needed to derive load models that simulate different human movements as sole individuals or as individuals in a group. This research takes a fundamental look at the problem of modeling individual forcing functions and characterizing the correlation between individual forcing functions. 3.0 EXPERIIVIENTAL MEASURENIENTS 3.1 Experimental Devices The goal of the experimental measurements was to accurately measure dynamic force components due to specified human activities. Data taken from tests performed on a force platform and a floor system were obtained from three sets of experiments, two Figure 3.1 View of the Force Platform, Saul et. al. (1990) designed to measure forces due to the action of an individual; and the third designed to measure the force due to the action of a group. Two different test devices were utilized. One was a 4 ft by 8 ft force platform designed and instrumented as described by McDonald (1984) and Ebrahimpour and Sack (1988) used for measuring the vertical and 18 19 horizontal components of loads produced by individuals and small groups up to five participants, Figure 3.1. The platform has an approximate fundamental natural frequency Figure 3.2 The Floor System Lab, Saul et.al. (1990) of 20 Hz. which provides a flat transfer function so that imposed loads are transmitted to the sensors without distortion in the frequency range of the load; this is essentially the definition of a force platform. The platform is supported vertically by five transducer assemblies (load cells), and horizontally by three in each direction. The other device, as shown in Figure 3.2, is a 12 ft by 15 ft floor system designed to measure total vertical loads generated by a group of people, up to forty participants, Ebrahimpour (1989). The floor system has a fundamental frequency of about 26.4 Hz and was instrumented with load cells and linear variable differential transformers (LVDTs) to measure the vertical component of human loads,Vankleek (1988). More detailed information on the data acquisition hardware is provided by Sack and Ebrahimpour (1988), and Burke et. a1. (1985). 20 Data-Set N0. 1 [101 [85] Transient Periodic dom Single Jump Penodc Jumping Random Jumping [26] (431 [121 L l i 1-Person 2-Persons 5—Persons 1-Eerson 2-Persons [10] l§L El [5] 1 arson 2-P ersons anions [$01 1201 [91 l l T T l 7 1 Hz 2 Hz 3 Hz 5 Hz 1 Hz 2 Hz 3 Hz 5 Hz [5] [51 1'5) 15] [21 (2] [21 171 m. 2‘... l... in. [51 [51 [5] l5] Data-Set N 0. 2 13511 J Transient Random Jumping PerLdle I154: I391 [“181 L Single Jrimp Sitting down Standirjg up 114:: 2-Per.r(s) 4901(5) [58] 1481 1481 l 191 15] 1-Per. 2-Per. (s) 4-Per.(s) 14:1“ 2433“! s) 44394, 5) 1451 I91 [41 I351 191 I41 1 , —Per. 5) 4-Per.(s) [351 [4] Jam rrig Jouncing Shalying [1 81911“! (“I41 F I l T It 4H1l l l 2 Hz 3 Hz 4H: 3 Hz 4 Hz [5°] [49] [49] 1-Per. MZ-PerSs). 4—P‘esr)(s) hiflZ-Pl 5) 44°19?) l——T—_——‘|l‘———1—-—l 1-Pelr Mtg} 4—Per(s) 1-gerj 2-Per(s) 4-Per(s).1-Per2-Per(s).4-Per(s).1-Pler. 2—Per(s). 4vPer(s) 1-Pler. 2-Per(s)4-Pen(ls) r481 :48: [49] 931.243 s). 4-P Ws)|p£—1——j1-PerJZ-Pegsl Mafia) Data-Set N0. 3 [11"] V V 1 sale are warm 1 I i l T 14’ . 9-P . 10-P . 20-P . 30«P . [546]! [9? (s) [631(5) [1 8elds) [g]! (S) Figure 3.3 Human Load Experimental Data-base 40-Per.(s) 151 21 In contrast to the force platform which by definition is designed so that input equals output, the floor system is modeled as a nine-degree of freedom vertical dynamic system where the output is read as vertical displacement of the floor. 3.2 Data Collection Dynamic human loading can be categorized as transient, periodic, and random. Transient loads are temporary non-repetitive loads from which the structural response is damped out before the next transient action is expected to occur, McDonald (1984). Types of transient actions could include a single jump, suddenly standing up, or falling into a seat. Human load modeling in this study is limited to the transient and periodic motions. 1 45 As illustrated in Figure 3.3, large Single Jump 2 19 74 number of tests, were conducted on both 4 4 5 6 the force platform and the floor system for 1 35 individuals and groups of up to 40 Sitting down 2 9 48 participants. Prerecorded pulses (prompts) 4 4 were played at the desired rate through Standing up 1 35 48 2 9 loud speakers, and the participants were 4 4 requested to perform a specific action at the pulse rate. The output from these tests were utilized to derive and verify the individual and group effect of human movements. Table 3.1 summarizes the collected 22 transient actions by the type of action (single jump, sitting down, and standing up). In Table 3.2, the periodic actions are categorized by the stimulus frequency for the motions in the study (jumping, jouncing, and swaying). Table 3.2 Periodic Actions Data-base 23 LSelect motion I< L Iread parameters I L +I New File I l I Read Raw Data ] I I Find Force-Historiesj N0 Motion is No Periodic Transien Random I Find Peaks I I Find Max. Peak I J r { Compute Statistics I‘ l I Save forces I We w L Plot Forces I L Calculate Autocorrelation and Autospectra I Plot Spectra I T —>I Save to Output I -———>J No (Save Statisticsj I Calculate Average Spectral Data ast No Motion Figure 3.4 DATA.FOR Flow Chart 24 3.3 Data Processing The raw data was collected in digital format as voltage readings. To use this data for the development of load models, a FORTRAN program, "DATA.FOR" was developed. For each motion type, the program reduces the transducer data into force— time series in the three principal axis of motion; vertical, and two perpendicular horizontal components Then the program calculates the statistical parameters of each force history, plots the forces in the three directions, computes and plots autocorrelation and autospectra for periodic motions, and calculates the average autocorrelation and average autospectra for force-time histories considered in the analysis. A listing of the program is included in Appendix A. Figure 3.4 is a flow chart to illustrate the components of the program "DATA.FOR". Input to the data processing program consists of three sets of files; a) raw data files, b) processing files, and c) calibration files. The raw files contain the voltage readings and information about the type of test, weight of person or persons performing the experiment, and time and date of the test. The processing files include "PARAM.DAT" which is a collection of parameters needed to process the raw data and to compute the statistical and spectral information. Font files which were needed to plot graphs on the monitor are also a part of the processing files set. Calibration files hold the calibration factors needed to convert voltage into forces. These calibration factors were computed in an earlier stage of the research project and the process of obtaining these factors are documented in Ebrahimpour (1988), and VanKleek (1989). 25 The program saves the computed statistical information and spectral analyses in separate files which have the same file names as the raw data, but with different l SJT1A010.FRC Force Forces in the x1, and 2 directions 2 SJT 1A010.ACN Autocorrelation Autocorrelation for XJ, and 2 directions 3 SJT1A010.ASP Autospectra Autospectra for the three directions 4 SJ-SUMM.OUT Final output Statistical information 5 SJ—AVG.ACN Correlation When more than one file grocessed 6 SJ—AVGASP Spectra When more than one file processed extensions. For example, a raw data file name "SJT1A010.RAW" which contains data from a single jump (SJ) performed as a transient action by one person (Tl), from data set (A), and was saved as test number (010) of the raw data collection (.RAW) would produce the files shown in Table 3.3 when processed by the program "DATA.FOR": Output to the monitor is composed of force-time plots and spectral analysis plots (for periodic actions). Figures 3.5 and 3.6 are views of typical output from DATA.FOR. Maximum peak values in the three perpendicular axis and the weights of the participants are among the statistical information reported on the figures. The vertical forces measured using the floor system were computed from transducers voltages using the Fortran program "FLOOR.FOR". The listing of the program is provided in Appendix A. The program provided an effective and accurate method to process the test data gathered from the floor system. The procedure employed in the program to convert LVDTs voltage readings to vertical forces is explained in detail by 26 DATAF O R Output FXAS A name to WEIGHT w as A RATIO ro wsraur .30 .08 r: .25 F .04 O O a 2 3° 2 a E .15. E .00 . n '10 n on A ~05. A. ' .7 .00' '7 '-°‘ 0 .35 0 “£3 -.10 -.09 as. -.10. FILE I PJ318002.RAW WITH 1 PER. F g PERIODIC MOTION - PJ TEST - 002 E l Q ; weights - 140. 0. 0. 0. 0. n ,5 g g Total-- 140. ES'rrrtxrz- 143.17 lbs A ‘ 'i i l i 1T ' ' m ' 3’ L” " 6 or 15 max ram: 0 " LUII L X-dir .316 3.912 Y-dir -.090 3.610 30 do Z-dir 3.134 3.618 TlMEUN-I 211' marooottrnru.... Figure 3.5 Force Type History (DATA.FOR Screen Output) DATA. F O R Output SPECTRAL ANALYSIS AUTOCORRELATION AUTOSPECTRUM 1.00 1.29 _ .eo« 1.00« f: a t c .901 3 -40‘ 5 o . 1 .201 ' .601 f 0 3 .00" 3 i s .40« l :3. -.20- ; .zoi ’-“°I .o 20.0 :0? 303 33.070110 120.0 .000 rfo in 3:0 {.0 6.0 do 7:0 $5.0 3:01:10 Les Fmfl-ltl FILE: PJ313002 HMS: 1.1612 FEAK: 3.18 MEAN: .06 Figure 3.6 Spectral Analysis (DATA.FOR Screen Output) 27 VanKleek, ( 1988). The floor system was modeled as a nine degree-of-freedom system and the program FLOOR.FOR was used to solve the following equation of motion: mm [M] [A] + [61m + [K] [x] -------- [3.11 where : FM : applied nodal force vector M9,, : mass matrix AM : acceleration vector C9,, : damping matrix V9,, : velocity vector K M : stiffness matrix 179,, : displacement vector LVDTs Readings were converted from voltages to displacements using the proper calibration factors. Displacement histories were differentiated numerically to obtain velocity and acceleration. The vectors were used in Eq. [3. 1] to calculate the force vector 17,. The total force-time history on the floor were computed by summing the force vector at every time t. Input to FLOOR.FOR consisted of the raw data files, calibration matrix, and the predetermined mass, damping, and stiffness matrices. Output files included force-time history generated at each node location (LVDT location) on the floor, total force-time history, and a summary of statistics of force peaks. Figure 3.7 shows a typical plot of force-time output from FlOOR.FOR 3.4 Statistical Output Summary: The force history, or shape of the force function with time, for each of the actions measured, proved to be consistent. Since each action appeared to result in a reproducible 28 function, it was hypothesized that each function could be accurately modeled. Although each force history for the same action is similar in shape, each is statistically different; that is, the force functions vary within narrow bounds. As a result an accurate model would be probabilistic. Statistical Results from the data processing program DATA.FOR output are stored on magnetic media for later use in the analysis. They are arranged by motion type (Transient, Periodic, and Random). For each motion type, output forces are normalized with respect to the participants' weights and referred to as ”Dynamic Force Ratio ”. Absolute maximum dynamic force ratios and actual and estimated static weight are reported for an individual and for a group. Following are abbreviations that apply to values in these tables: Fxmax : absolute maximum force ratio (to individual's weight) in the x- direction. Var : sample variance STD : sample standard deviation C.O.V. : coefficient of variation Table 3.4, a typical statistical table, shows the absolute maximum dynamic force ratios by one individual performing a single jump in-placc. It can be seen from the table that force ratios in the vertical direction are significantly more than the other two components in the horizontal direction. The table shows also that an individual would generate an average vertical force of about 3.19 times his weight when performing a single jump. The minimum vertical dynamic force ratio observed in this collection of tests was 1.54. Thirty five males and fourteen females participated in this experiment. 29 Their average weight was 161.36 lbs. with a standard deviation of 30.80 lbs. and weight range from 103 to 242 lbs. 251-342;, if,- _ . 7 30 For each motion type, a table summarizing the statistical information of individual and group force ratios is prepared for later analysis. Table 3.5 , for single jump motion, Table 3.5 Maximum Force Ratios (Single Jump) Fymax l 2 4 .60 .37 .24 .28 .21 .10 .21 .03 .04 .09 1.24 0.65 1.49 0.73 .31 .11 0.03 . 0.06 0.03 0.01 0.33 .17 . 0.25 .18 0.09 45.79 102.3 63.78 .41 10 Persons (Perlodlc Jumplng at 2 Hz.) ll /\ fl «Li ill -- f \1 U \/ if \f U \«x’ W V 1.5 2 2.5 0.5 Dynamic Force R0110 o I d 0 0.5 Time (see) Figure 3.7 A Typical Force-Time History (Floor.For Output) presents a typical sample of such tables. It shows that an increase in the number of people performing the same action at the same prompt tends to decrease the overall dynamic force. A single jump exemplifies this phenomena since it is difficult for typical individuals to maintain coherency in a pulse action like a single jump. 31 Tests performed on the floor system were analyzed for peak statistics using FLOOR.FOR. Some of the statistics are tabulated in Table 3.6. Data collected from the floor system are limited in two aspects; only the vertical component of the dynamic force was recorded; and data collection is limited to only periodic jumping at 2Hz. It was observed that individuals could maintain better coherency at this frequency range. Data collected from the floor system were used in this study to verify a load model hypothesized to simulate a large crowd. Table 3.7 shows that the larger the group that performs the same action at the same time, the harder it is for individuals to maintain coherency. The maximum dynamic force ratio decreased from 4.54 for one persons to 1.14 for forty persons performing periodic jumping at 2 Hz. Figure 3.8 shows the group dynamic ratios for groups of individuals up to forty. The sample size of the four participants was the smallest among all other groups. Therefore, the standard error of the mean for that particular group shows a 32 wider range than all other groups. The figure also suggests that the relationship between the number of participants and the dynamic force ratios may not be linear. Table 3.7 2 Hz.-Periodic Jumping Peak Analysis Group Dynmalc Force Ratios (Periodic Jumping, 2 Hz.) 5 .00 4.50 4.00 Dynamic Force Ratio !° N S" 9’ 0 or 0 or o o o o -n 01 o 1 .00 0.50 0.00 1 2 4 1 O 20 30 40 No. of Participants Figure 3.8 Mean, Min, Max Peaks For 2 Hz. Periodic Jumping Since periodic jumping at 2 Hz. was found to produce the maximum coherence among participants, it would constitute an upper limit of dynamic load ratios of one 33 individual in group. A non-linear regression model is proposed to estimate "Dynamic Group Factor " DF; an empirical factor to generate dynamic force ratios as function of group size n . The model is in the form : where regression coefficients b, m are estimated by transforming Eq. [3.2] to a linear form and solve the linear regression, Chapra (1988): Log(DF)==nLog(m)+I.og(b) .......... [3.3] D F = b * m n .................. [3.2] The regression model was solved to estimate the regression coefficient for mean, minimum, and maximum dynamic factors. Figures 3.9 presents the regression model in comparison to the measured values. The figure also shows that the range between the Dynamic Factor Ratio (DP) (Exponential Mode 4 i» * "m" min 3 5 IN / max .I AWW [fathered Maxinurn Peekd + actual mean 3 _ 1‘ W \ —I-— actual min 2:: 2 5 3‘ \ m A aotualmax a ' ' ‘N 33 N .9 \ s 2 -- s 3‘ I 1.5 ir‘ LIIL llmlll—rlll—l—ilWlWlllililT‘lTilimllWTllillllliTl—Tlllflfliliil—Fllll—WTIIlilTfllllllITlll—T‘l 0 10 20 30 4O 50 60 7O 80 90 NeolWln) Figure 3.9 Log-Linear Regression Model 34 maximum and minimum dynamic group factors decreases as the crowd size increases. This formulation could be used as a tool to estimate an upper bound for the mean and range of the dynamic force ratios of any group size n. The pmposed Log-linear Single Jump A single jump in place from a standing position I Standing Up Standing up from a sitting position ]I I Sitting Down Dropping into a seat from a standing position I Jumping Periodic jumping in place I Jouncing moving up and down without leaving the floor ..... LSWaying Swaying side to side from a sitting position regression can be simplified to take the form: DF = a * B" .................. [3.4] where : or is 2.98, 1.73, and 4.08 for the mean, minimum, and maximum Dynamic Group factors and B is a constant factor equals to 0.97. Although several other regression forms were examined, the proposed regression in Eq. [3.4] proved to be the best—fit model. Nevertheless, it is important to emphasize the preliminary nature of the proposed regression. It is beyond the scope of this research to further investigate the theoretical aspects of the regression models in spite of their practical importance. Further statistical analysis could provide a better and more representative form to describe DF. 4.0 MODEL ANALYSIS 4.1 General In this part of the research, the focus is analyzing human loadings data to include and quantify their random nature. For each transient action in the study, a method based on passing straight lines through control points that define the force-time history is used to simulate human loading. Best-fit distributions for the control points were determined. Periodic loadings were analyzed as a series of impulses where the forcing function is divided into cycles, each of which is a full period of the motion considered. Statistics of each impulse in the force-time history were computed and analyzed for correlation and peak distribution. 4.2 Methodology The dynamic force resulting from a single activity, such as a jump, is a time varying function F,(t) graphically represented as a single pulse with force changing with time over a relatively short interval, Fig. 4.1.a. This dynamic load will vary somewhat from individual to individual and from one time to the next for the same individual. In addition, when two or more individuals attempt the same action at the same time, the resulting loads will vary somewhat from one another but may be influenced by each 35 36 other, Figure 4.1.b. The function F,(t) due to the action of a single person i may be a pulse or a continuing response. Table 4.1 summarizes motions considered in the study: Table 4.1 Human Movements Considered in The Study Single Jump A single jump in place from a standing position Standing Up Standing up from a sitting position II . I Sitting Down Dropping into a seat from a standing position I Jumping Periodic jumping in place [Jouncing moving up and down without leaving the floor Swaying Swaying side to side from a sitting position For each movement type, the forcing function is assumed to be a function of several independent random variables such as time, response lag, and group effect. The parameters of the function due to an individual performing a predefined human movement can be represented as F,(t): 5,0), = w l 1 +fz (t, ¢,f)l ............. [4.1] Fin), = w j, ( t, o, f) ............... [4.2] F,(t)y = w fy ( t, ¢,f ) ............... [4.3] where the random variables are: w : Weight of individual t : Time in seconds (I) : Response Lag (Time lag between the prompt and the individual reaction) f : Frequency of motion (not a factor in transient actions) 37 The dynamic components of the forces due to one person is approximated by passing line segments through control points defined as shown in Figure 4.2 for the case of a single jump, Tuan (1985). Periodic motion is modeled as a series of impulses, where the shape of each impulse is defined by passing line segments through control points as shown in Figure 4.3. The mean and the variance of arrival time and amplitude is determined for each control point. The probability density functions for the best-fit distribution is evaluated for arrival times and amplitudes of control points. To simplify the modeling process, several assumptions were made. They are: 1.0 0.6 ~ 0.2 I -0.2 '1 -05 l Dynamic Ratio .00 1.00 2.00 3.00 4.00 ‘ .00 1.00 2.00 SOT 4.00 Tlmslsn) Time(sec.) a) One Individual 0) Three Individuals Figure 4.1 A Single Jump l Initial computations are confined to the vertical component of the dynamic force. 38 A Single Jump, One Person A J‘A‘ —P V-“ A w m m m m o 2 2 1 1 73:88“. Time ( sec. ) Time ( sec.) 10 Persons Periodic Jumping Figure 4.2 Control Points for A Single Jump Figure 4.3 Control Points for a Periodic Motion -50 4 ~100 a -150 39 I Random processes used in the analyses were assumed to be adequately characterized by second-order models (models that can be sufficiently defined with the mean and the standard deviation). I To account for participants weight variabilities, a general model based on the statistics of the weight distribution of the US . p0pulation were assumed for the study. I For periodic motions, A model assuming that the person is always trying to maintain coherency with the prompt was assumed. 4.3 Transient Loading The study includes single jumps, standing up from a sitting position, and sitting down. A typical transient forcing function F,(t) has the following random variables: Response Lag (I) : A random variable that represents the time between the prompt and the response of an individual. Available data were utilized to find the probability distribution function of the Response lag, and to determine whether it depends on the type of the transient motion. Amplitudes a 1,, a,,,..., am: Amplitudes of the control points that best—represent the transient function considered were determined. Covariance matrix and correlation with individual's weight were computed and analyzed. Arrival Times tfi, thy”, t": Similar to the amplitude coefficients, statistics on arrival times were computed. The first arrival time t,,. for an individual i, is always assumed to be equal to the response lag (I) assumed for the individual. 40 4.3.1 Single Jump The single jump motion is used to explain in details the process of analyzing the measured data. Other types of transient actions considered in this study (Standing up, and Sitting down) are analyzed similarly. Differences in the analyses process are explained later in this Chapter. For each measured force—time history, as a stochastic process, the weight of the person performing the test was removed from the process by subtracting it from all Vertical Force Ratio (Single Jump, One Individual) .2 E i ll. 0 'E i o 0.5 1 1.5 2 2.5 Time (See) Figure 4.4 A Typical Single Jump force-time series points. The process was normalized by dividing each reading by the weight of the participant. Figure 4.4 shows a typical vertical component of a single jump. 41 Eight control points were assumed sufficient to describe the vertical component of the forcing function, as shown in Figure 4.4. The prompt marks the beginning of the process. The first control point denotes the start of the individual's motion. The time lag from the beginning of the process to the first control point is defined as the Response Lag 0. The individual is in a static position at this portion. The process consists of three segments; first is the preparation for the jump, from control point 1 to 4; then airborne time between control point 4 and 5; and finally the landing portion where control point 6 marks the peak landing. The process ends at control point 8 when the person returns to the static position, Figure 4.4. A program, CONTROL.FOR was developed to automate the process of analyzing the measured data. Input to the program consists of a parameter file and the collection samples of force-time histories for the chosen action. CONTROL.F OR proceeds as follows. I Reduce the force-time history to 2.5 seconds after the prompt. It has been observed that 2.5 seconds is well beyond the duration of any action period considered herein. Force-time histories were aligned such that the prompt marks the beginning of the file. I Differentiate the reduced function with respect to time to compute the slope function. The slope function is used to find the control points. I Find control points using two criteria. First, the change in the slope direction will define peak control points 2,3,6, and 7. Control points 1,4,5, and 8 are defined as points where the slope values changes significantly. The program compares the 42 slope before and after each point and determines whether it exceeds a preset range value. The range value establishes the second criteria. The response lag (I) is determined as the time from the prompt to the first control point. The response lag is then removed from the process and treated as an independent random variable. This assumption is verified later in the analysis. Arrival times are then adjusted so that actual dynamic motion starts at time g=aa Reduced force-time history data along with amplitudes and arrival times of the control points are saved in a file which has the same name as the input force-time history with an extension ".CRL", The Subroutine PLOTXY, is then called to plot the reduced function and control points on the same axes to visually ensure that control points were properly captured. Arrival times and amplitudes of control points are also shown on the screen. The program repeats the previous steps for all samples of force-time histories in the selected human motion (a single jump in this case). An output file (SJ -CRL.OUT) saves arrival times, amplitudes of control points, and peak values in the three principal axes from each sample in the collection of files processed by the program, Table 4.2. Subroutine MYSTAT is called to compute the important statistics of control points amplitudes and arrival times, Table 4.3. The computed statistics are saved in the output file. 43 Table 4.2 Control Points for Single Jump WT (I) T1 T3 T4 T5 T6 T7 T8 A1 A2 A3 A4 A5 A6 A7 A8 M“ 165 0.18 0.12 0.32 0.44 0.85 1.00 1.29 1.62 -0.10 -0.64 2.17 -1.00 -l.00 1.41 -0.45 0.00 172 0.15 0.15 0.41 0.56 0.97 1.03 1.41 1.53 -0.04 -0.68 1.92 -0.97 -0.98 2.55 -0.57 0.09 162 0.21 0.09 0.41 0.59 1.06 1.15 1.79 2.00 -0.18 -0.54 1.54 -l.03 -l.01 1.44 -0.45 0.02 181 0.29 0.09 0.35 0.44 0.88 0.97 1.27 1.38 -0.15 -0.45 2.24 «0.95 -0.92 2.38 -0.51 0.21 175 0.15 0.00 0.21 0.35 0.82 0.94 1.38 1.62 -0.05 -0.05 2.07 -0.95 -0.97 4.24 -0.20 0.03 165 0.21 0.12 0.27 0.38 0.74 0.79 1.18 1.35 -0.06 -0.74 2.42 -1.01 -l.00 2.69 -0.53 0.02 172 0.15 0.12 0.32 0.44 0.85 0.91 1.24 1.59 -0.04 -0.98 2.25 -1.01 -l.01 4.36 -0.49 -0.03 162 0.18 0.15 0.41 0.68 1.15 1.24 1.91 2.09 -0.08 0.52 1.29 -1.00 -l.01 2.04 -0.51 0.10 181 0.32 0.09 0.32 0.44 0.79 0.88 1.09 1.24 -0.07 -O.62 1.86 -l.00 -l.01 3.24 -O.7l 0.07 176 0.27 0.00 0.09 0.18 0.53 0.59 0.85 1.53 -0.03 -0.03 2.56 -l.02 -l.07 3.86 -0.79 -0.25 196 0.27 0.09 0.47 0.59 0.97 1.06 1.29 1.47 -0.03 -0.47 1.61 -0.93 -0.97 3.49 -0.34 0.20 140 0.50 0.06 0.35 0.50 0.94 1.06 1.24 1.27 -0.23 -0.41 1.70 -0.92 ~0.94 2.62 -0.11 0.03 186 0.53 0.09 0.65 0.77 1.32 1.38 1.82 1.94 -0.06 -0.38 1.41 -l.03 ~0.97 3.43 -0.43 0.11 150 0.32 0.09 0.24 0.38 0.77 0.85 1.12 1.24 -0.08 -0.72 2.65 -0.89 -0.95 4.60 -0.32 0.03 189 0.35 0.12 0.38 0.53 1.00 1.09 1.27 1.44 ~0.07 -0.65 2.06 -0.98 -0.99 5.04 -0.74 0.07 163 0.24 0.18 0.47 0.65 1.09 1.18 1.50 1.62 -0.01 -0.56 1.60 -0.91 -0.92 1.95 -0.45 -0.14 125 0.38 0.09 0.38 0.53 0.94 1.00 1.47 1.68 -0.07 -0.49 1.56 -0.89 -0.90 5.04 -0.18 0.03 234 0.35 0.09 0.29 0.47 0.65 0.77 0.97 1.06 -0.06 -0.29 1.46 -0.85 -0.88 2.12 -0.23 0.06 132 0.29 0.12 0.27 0.38 0.85 0.94 1.21 1.50 0.02 -0.72 2.93 -0.88 o0.87 5.94 -0.58 -0.12 130 0.24 0.12 0.27 0.41 0.71 0.82 1.00 1.18 -0.04 -0.60 2.13 -0.95 -0.94 2.43 -0.34 0.00 146 0.27 0.12 0.24 0.35 0.71 0.82 1.00 1.18 -0.01 -0.85 2.47 -0.90 -0.96 2.31 -0.28 0.19 242 0.18 0.12 0.44 0.56 1.03 1.09 1.56 1.71 -0.07 -0.43 1.31 0.98 -0.98 3.49 -0.42 0.01 184 0.35 0.12 0.38 0.50 0.94 1.03 1.27 1.47 -0.06 -0.31 2.44 -0.98 -0.97 4.05 -0.53 -0.03 155 0.35 0.18 0.44 0.59 1.09 1.15 1.44 1.62 -0.05 -0.41 1.84 -0.96 -0.96 3.12 -0.69 0.06 200 0.21 0.18 0.47 0.59 1.00 1.09 1.29 1.50 -0.01 -0.36 1.86 -0.97 -l.02 4.47 -0.36 -0.04 191 0.21 0.09 0.21 0.35 0.56 0.62 0.82 1.03 -0.03 -0.70 1.86 -0.93 -0.93 3.82 -0.54 -0.09 166 0.35 0.09 0.27 0.41 0.71 0.79 1.06 1.18 -0.11 -0.49 2.08 -0.90 -0.95 2.13 -0.56 0.08 215 0.38 0.09 0.27 0.38 0.79 0.85 1.12 1.24 ~0.01 -0.54 2.27 -0.90 -0.87 2.84 -0.63 0.13 154 0.32 0.12 0.27 0.41 0.77 0.85 1.15 1.29 -0.04 -0.82 2.20 -0.98 -0.97 4.24 -0.53 0.07 153 0.21 0.12 0.35 0.50 0.97 1.06 1.29 1.44 -0.03 .0.62 1.99 «0.90 -0.90 2.97 -0.63 0.43 164 0.18 0.06 0.27 0.38 0.79 0.85 1.12 1.59 -0.10 —0.59 2.79 -1.01 -l.01 8.85 -0.40 0.01 120 0.29 0.09 0.24 0.38 0.56 0.62 0.77 0.79 -0.05 -0.42 1.29 -0.92 -0.99 1.85 -0.13 0.04 103 0.27 0.09 0.21 0.32 0.50 0.59 0.77 0.79 -0.05 -0.72 2.23 -0.88 -0.90 2.80 -0.15 0.01 126 0.27 0.00 0.06 0.18 0.32 0.47 0.74 0.85 -0.04 -0.04 1.34 -l.01 -l.10 1.56 -0.53 0.01 116 0.21 0.18 0.38 0.50 0.88 0.91 1.21 1.38 0.04 -0.41 1.94 -0.85 -0.86 1.77 -0.59 0.49 185 0.24 0.12 0.38 0.47 0.85 0.91 1.21 1.56 0.05 -0.38 2.07 -0.94 -0.96 4.27 -0.34 0.02 145 0.35 0.09 0.38 0.53 0.88 0.97 1.32 1.38 -0.11 -0.43 1.38 -0.82 -0.93 2.03 -0.08 0.01 140 0.24 0.06 0.18 0.27 0.50 0.59 0.85 1.00 -0.07 -0.57 2.57 -0.69 -0.86 2.23 90.28 0.04 170 0.29 0.09 0.32 0.47 0.82 0.91 1.12 1.24 0.01 -0.42 1.60 -0.87 -0.89 3.06 -0.46 0.01 120 0.41 0.06 0.21 0.29 0.65 0.74 0.94 1.12 -0.22 -0.58 2.70 ‘0.83 -0.93 3.54 -0.69 0.21 171 0.27 0.09 0.27 0.38 0.68 0.77 1.06 1.35 0.03 -0.78 1.66 -0.86 -0.87 1.95 -0.32 0.20 177 0.24 0.12 0.27 0.41 0.65 0.74 0.94 1.09 -0.03 -0.49 1.65 -0.87 -0.91 2.72 -0.47 0.14 124 0.29 0.06 0.32 0.47 0.77 0.88 1.12 1.06 -0.21 -0.59 1.32 -0.93 -0.96 2.13 -0.37 -0.29 128 0.32 0.06 0.38 0.53 0.82 0.91 1.12 1.38 -0.16 -0.46 1.06 ~0.79 -0.86 1.89 -0.05 0.02 110 0.24 0.12 0.29 0.41 0.71 0.79 1.06 1.24 -0.05 -0.41 1.71 -0.89 -0.94 3.04 -0.32 0.02 44 I All output files are closed and the program report the number of processed files and any errors during the run. Output from the program is used next to study the correlation among control points and to find the best-fit distribution for arrival times and amplitudes of each control point. Correlation Analysis In order to better understand the relationship between the control points arrival times and amplitude ratios, a correlation analysis was conducted. The correlation coefficient r was calculated for all control points pairs and the results are tabulated in Table 4.4. The coefficient of correlation is a measure of the degree of linearity in the relationship between any two variables X, Y, Pollard, (1979). The correlation coefficient is calculated from the covariance matrix as follows: N COV(x,y) =—;i£(xi-xm)(yi—ym) ......... [4.4] r = ———"°"("’y) ................. [4.5] x” oxoy where cov (x , y) is the covariance matrix between variable X and variable Y, N is the number of samples, x, is the value of the X variable from sample 1', and x, and o, are the mean and the standard deviation of X. The coefficient of correlation takes values from -1 to +1 where the negative values suggest negative correlation (high values of x are associated with low values of y) and positive values indicate positive correlation (high values of x are associated with high values of y). The perfect positive linear correlation will result in a value of 1 for the coefficient. An r value of zero does not imply 45 independence between .1: and y. It only indicates lack of any linear relationship between x and y. In general, 1001’ gives the percentage of the total variation of the variable Y which is explained by, or is due to, the relationship to variable X, Freund (1980). It is however important to emphasize that r measures only the strength of the linear relationship and that high correlation values does not necessarily imply a cause-effect relationship. If the considered variables were drawn from a bivariate normal distribution, the correlation 0.10 0.45 0.82 0.90 1.19 1.37 -0.0 -0.5 1.93 -0.93 -0.9 3.16 -0.43 0.05 103 0.15 0.00 0.06 0.18 0.32 0.47 0.74 0.79 ~0.2 -0.9 1.06 -1.03 -1.1 1.41 -0.79 -0.29 242 0.53 0.18 0.65 0.77 1.32 1.38 1.91 2.09 0.05 -0.0 2.93 -0.69 -0.8 8.85 -0.05 0.49 139 0.38 0.18 0.59 0.59 1.00 0.91 1.18 1.29 0.28 0.95 1.88 0.34 0.24 7.45 0.74 0.78 949 0.01 0.00 0.01 0.01 0.04 0.04 0.07 0.08 0.00 0.04 0.21 0.01 0.00 1.88 0.03 0.02 31 0.09 0.04 0.11 0.12 0.19 0.19 0.26 0.29 0.06 0.20 0.46 0.07 0.06 1.37 0.18 0.13 . 19 30.96 41.47 33.22 25.91 23.41 20.70 22.03 20.82 -100 -37. 23.74 -7.50 -5.7 43.4 -42.5 264.2 @3555; 5 0.01 0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.01 0.03 0.07 0.01 0.01 0.21 0.03 0.02 coefficient can be further examined by using the Fisher Z transformation which estimates the confidence intervals for calculated coefficient of correlation. The test is based on a change of scale from r to Z which is given by: where the distribution of Z is approximately normal and a direct function of the true population coefficient of correlation p. The confidence interval for p is constructed from known statistical tables and the following (1-01) confidence interval: 46 Z-_zLfl— .. « ...-- :: i ; z: . .. @1082; 158. 5714 34. 7261 0.0922 0.0895 LOG-NOR 0.2521 0.1088 0.1672 0.3588 NORMAL 0.2445 0.1035 0.1448 0.087 LOG-NOR 0.4631 0.1372 0.1354 0.0895 LOG-NOR 0.5328 0.1523 0.1737 0.128 LOG-NOR 0.6025 0.1568 0.217 0.1614 LOG-NOR 0.7479 0.1713 0.146 0.1001 LOG-NOR 0.8622 0.1746 0.158 0.1272 LOG-NOR 0.0417 0.0544 0.2608 0.2315 LOG-NOR -0.7595 0.6027 0.2467 0.1394 NORMAL 0.8298 1.3254 0.3551 0.0847 LOG-NOR 0.8609 1.3101 0.342 0.1299 LOG-NOR 2.1178 1.7622 0.224 0.0905 LOG-NOR -0.3263 0.2298 0.2434 0.1171 NORMAL 0.1072 0.0816 0.1857 0.0817 LOG-NOR 0.2497 0.1541 0.1555 0.1033 LOG-NOR 0.2608 0.184 0.2095 0.0899 LOG-NOR 0.0766 0.0738 0.3124 0.1689 LOG-NOR 0.0808 0.0672 0.2444 0.164 LOG-NOR 2.1 178 1.7622 0.224 0.0905 LOGoNOR 4.3.3 Standing up The vertical component of the forcing function for the standing up motion was found to be sufficiently approximated by six control points, Figure 4.7. The program CONTROLFOR was modified to analyze the motion. Statistics for control points and peak values are shown in Table 4.9. The table shows that the average peak force in the Table 4.9 Standing Up Control Points Statistics 54 Vertical Force Ratio (Standing up, One Individual) Dynamic Force Ratio 0 0.5 1 1.5 2 2.5 Time (See) Figure 4.7 A Typical Standing up .......................................................................................... ........................................................................................................................................ 158.5714 34. 261 0.0922 0.0895 0.3135 0.1576 0.1319 0.1196 LOG-NOR 0.2344 0.1113 0.115 0.1298 NORMAL 0.5285 0.1232 0.0985 0.0884 LOG-NOR 0.621 0.1568 0.1199 0.1335 NORMAL 0.7941 0.2276 0.1621 0.1096 LOG-NOR 0.9471 0.2416 0.1205 0.0799 LOG-NOR 0.0159 0.0088 0.1104 0.1623 NORMAL 1.0104 0.8999 0.2032 0.1074 LOG-NOR 0.9604 0.8472 0.308 0.1619 LOG-NOR 0.9419 0.846 0.305 0.1574 LOG-NOR 1.0809 1.2655 0.2624 0.1917 LOG-NOR 0.0805 0.0893 0.2452 0.1862 LOG-NOR 0.3708 0.3255 0.2508 0.136 LOG-NOR 0.5074 0.4085 0.3124 0.1742 LOG-NOR 0.1071 0.1523 0.3053 0.1304 LOG-NOR 0.145 0.1323 0.1946 0.1015 LOG-NOR 1.4718 1.2511 0.229 0.1477 LOG-NOR 55 vertical direction for a person standing up from a sitting position was about 1.47 times his or her weight. Force ratio components in the orthogonal horizontal directions were 0.5 for the average force perpendicular to the seating and 0.14 for the average parallel force ratio. Comparison with same statistics for the sitting down motion reveals that the later produced larger forces in the vertical direction. This observation is expected since the person standing up is moving against gravity while moving with it in sitting down. The same argument could also be used to explain the finding that the average time needed to perform the motion was 1.003 seconds with a delay after the prompt (response lag) of 0.291 seconds in comparison with only 0.86 seconds to perform the sitting down motion and a response lag of 0.25 second. Linear correlation coefficients were calculated for all control points pairs with the results tabulated in Table 4.10. Finding the best-fit probability distribution to represent each variable in the motion was accomplished by examining the K—S values shown in Table 4.11. The best-fit distribution for the response lag was found to be the Log-Normal distribution. This confirms the finding in Table 4.5 for the case of a single jump. The table reaffirms again that amplitudes can be best-represented by a Log-Normal Distribution. 4.4 Periodic Loading Periodic Loading includes: periodic jumping, jouncing, and swaying. The dynamic component of the vertical force due to one person performing a periodic motion is approximated by a series of impulses, where the shape of each impulse was defined by 56 Table 4.11 Correlation Matrix (Standing up) ............................................................................................................................................................................................................................................................................... '2'?" :3 . 0.339 -0.29 0.253 0.246 -0.00 -0.20 -0.33 0.359 -0.36 0.275 -0. 11 .017 .025 -0.19 .023 -020 -0.18 .015 -0.10 -0.17 -0.25 0.008 -0.38 0.341 0.344 -0.12 0.114 032 0.298 .032 0.158 -0.27 0.631 0.704 0.479 0.389 -0.28 0.441 -0.47 .048 0.285 -0.06 0.526 -0.46 0.488 .037 0.43 01.0736 0.629 0.566 -0.18 0.551 -052 -0.52 0.225 -0.10 0.388 -054 0.482 -027 0.40 .027 0.439 -042 -043 0.545 0.067 0.481 -042 0348 -0.60 0.59 .005 0.291 -019 .019 0.142 .005 0.243 -022 0.198 -0.31 0.25 5771-001 0.456 .0. 34 .035 0.501 -0.08 0.409 .037 0.248 .055 0.56 ' -0. 29 0.231 0.209 0.058 .010 -0.26 0.235 0.32 0.319 -009 ' -0.83 0.441 -0.01 0.804 .079 0.617 -0.36 0.70 ‘ . .044 0.113 -0.89 0.947 -0.77 0.429 -0.64 . 0.12 -0.89 0.949 -o.77 0.445 -0.66 “ -039 0.193 -0.62 0.92 -0.07 0.064 0.06 -0.17 -0.18 §§§§§§ .025 -015 -019 -0.10 -0.23 -0.17 40.20 .025 ;§; 0339 0.008 -0.28 -0.18 -0.27 -0.05 -029 -O.38 0.441 0.551 0.439 0.291 0.456 0.253 0.341 -047 -0.52 -042 .019 -034 0.246 0.344 -0.48 -0.52 -043 -0.19 .035 " 0.00 -0.12 0.285 0.225 0.545 0.142 0.501 . .0.20 0.114 -0.06 -0.10 0.067 -0.05 -0.08 -0.33 -032 0.526 0.388 0.481 0.243 0.409 -0.50 0.67 0.359 0.298 -0.46 -0.54 -0.42 -0.22 -037 0.537 .0. 59 .0. 40 0.45 -0.36 -0.32 0.488 0.482 0.348 0.198 0.248 0.275 0.158 -0.37 -0.27 -0.60 -0.31 -0.55 0.319 ~0.36 0.429 0.445 -0.62 0.064 -0.50 ' -0.11 -0.2710.431 0.406 0.592 0.254 0.566 -0.09 0.707 -0.64 -0.66 0.92 0.069 0.678 .0. 59 0457 . passing line segments through predetermined control points. Impulses were saved as series of transient actions with the following assumptions: I If the prompt frequency is within the physical limits of an individual to perform the prescribed motion, the individual in general attempts to synchronize the motion to follow it. It has been observed that when an individual is out of phase with the prompt in one cycle, an effort is made in the next cycle to adjust to it. Findings from the measured data asserted this observation. I When a person, due to physical limitations, is unable to maintain coherency with a high prompt rate (frequency higher than the range of 2-3 Hz.), Reiner (1987), The person would perform the action at the highest frequency rate possible with no attempt to adjust to the prompt. Although insufficient data were available to verify this assumption, visual observation of tests conducted at frequencies of 4 57 and 5 Hz. demonstrated the leaning of individuals to perform actions at their highest possible rate. I The dominant forces on the platform due to the side swaying motion is the horizontal longitudinal component (parallel to the seating). I It is assumed that individuals maintain the same level of activity throughout the length of the studied motion (approximately 10 seconds). I The studied motions are treated as weakly stationary processes. The validity of this assumption is discussed later in the chapter. I The analyzed motion is confined to the periodical part of the motion. Initial stages of the motion would be simulated separately as transient actions. This is necessary to preserve the assumption made earlier about the stationarity of the studied process. 4.4.1 Periodic Jumping Force-time histories were divided into cycles, where each cycle is a full period of the action considered. Three control points were sufficient to describe each cycle. One is the peak point in the cycle and the other two define the airborne time. Figure 4.8 presents a typical vertical force ratio of one individual jumping at 2 Hz. For each measured force-time history, the analysis was conducted using the program CONTROLFOR which was described in Section 4.3.1 for the case of single jump motion. The program was modified, however, to accomplish the following additional steps for analyzing periodic motions in this study: 58 : j. W..." l 1 1‘ l l J .0 0 J UUJU 0000‘ o 2 4 6 Time (Sec.) Figure 4.8 A Typical Periodic Jumping I The number of cycles was computed for each measured load-time history. The __ Dynamic Force Ratio starting time of each cycle was recorded along with the peak force in the cycle. I For each cycle, arrival time and amplitudes of control points were defined. Figure 4.9 shows the control points for the periodic motion in Figure 4.8. I An output file is created in which the periods of each cycle and control points data are saved. Table 4.12 is a sample of information reported in the output files. I Subroutine SPECTRAFOR is called to compute the autocorrelation, autospectrum fimction and the root mean square of the process considered. Typical screen output is shown in Figure 4.10 and 4.11 for the forces in vertical and parallel axes of motion. Output from CONTROLFOR shows that the average peak force ratio in the vertical direction for a person jumping periodically in place at a 2 Hz. prompted frequency was about 2.273 times the person 's weight, Table 4.13. This ratio decreased to 2.154 at 3 Hz. 59 Control Points Periodic Jumping (2 Hz., One Individual) '1‘ M 1111‘ ii 00/010111 1 111 0 11111111111111 J Time (sec.) Figure 4.9 Control Points for a Typical Periodic Jumping and 2.002 at 4 Hz. This finding shows that the higher the jumping frequency the lower the average peak force ratio a person generates. Maximum peaks also followed the same pattern and dropped from an average of 2.852 for 2 Hz. jumping to 2.750 and 2.448 for 3 and 4 Hz. respectively. At a prompt of 2 Hz. (0.5 sec. period), the period of the jumping ranged from 0.472 to 0.59 seconds with an overall average of 0.509 seconds. This means that on the average, individuals were jumping at a frequency of 1.95 Hz. while the prompt frequency was 2 Hz. This difference between prompt and response increased at 3 and 4 Hz. as shown in Table 4.14. The average response period for persons jumping at a prompt of 0.333 seconds (3 Hz.) was 0.36 seconds (2.77 Hz.) and for the case of a prompt period of 0.25 seconds (4 Hz.), it was 0.285 seconds (3.51 Hz.). This observation verifies the assumption made earlier about the inability of a person to maintain coherency with the 60 Table 4.12 2—Hz. Periodic Jumping 1 . 2 173 0.478 0.045 2.22 0.416 3.331 3 126 0.477 0.033 3.148 0.39 3.94 4 160 0.497 0.021 3.627 0.468 4.219 5 205 0.476 0 032 2 539 0.404 3.175 6 196 0.509 0.042 2.397 0.252 2.883 7 140 0.501 0.039 2.685 0.246 3.186 8 186 0.492 0.104 3.222 1.395 4.536 9 150 0.507 0.023 2.513 0.136 2.824 10 189 0 524 0.035 3 126 0.485 3.736 11 163 0 521 0.013 2 746 0.132 3.006 12 125 0 55 0.173 2 669 0.352 3.517 13 234 0 488 0 027 1 524 0.166 2.002 14 132 0 587 0 025 3 825 0.326 4.227 15 130 0 498 0 037 1 784 0.218 2.462 16 146 0 526 0 023 1 93 0.197 2.372 17 242 0 505 0 026 1 542 0.136 1.774 18 184 0 472 0 028 2 451 0.469 3.107 19 155 0 547 0 076 1 094 0.24 1.579 20 200 0 504 0 034 2 537 0.32 3.069 21 191 0 505 0 05 1 066 0.252 1.861 22 166 0 493 0 035 1 894 0.261 2.522 23 215 0 517 0 024 3 459 0.435 3.935 24 154 0 472 0 03 2 608 0.193 3.084 25 153 0 533 0 03 2 454 0.47 3.198 26 164 0.504 0 027 3 213 0.171 3.473 27 120 0.51 0.1 2 188 0.176 2.487 28 103 0 51 0.039 1 399 0.323 2.243 29 126 0 52 0.024 1 337 0.219 1.732 30 116 0 522 0 097 1 129 0.224 1.512 31 185 0 489 0 053 1 684 0.253 2.189 32 145 0 504 0 042 2 444 0.42 3.029 33 140 0 496 0 044 2 013 0.148 2.3 34 170 0 513 0.025 1 758 0.16 2.008 35 120 0 511 0.028 2 044 0.287 2.542 36 177 0 59 0.267 2 142 0.167 2.384 37 124 0 518 0 017 1 924 0.235 2.282 38 128 0 513 0 029 1 334 0.334 1.829 39 110 0 498 0 038 2 107 0.369 3.413 AVG 159.179 0.509 0.048 2.273 0.319 2.852 MIN 103 0.472 0.013 1.066 0.132 1.512 MAX 242 0.59 0.267 3.825 1.395 4.536 RNG 139 0.118 0.254 2.759 1.263 3.024 VAR 1135.89 0.001 0.002 0.498 0.044 0.651 STD 33.703 0.026 0.046 0.706 0.21 0.807 cov 21.173 5.113 95.406 31.053 65.763 28.278 SER 5.397 0.004 0.007 0.113 0.034 0.129 61 a .. Dos Ea PERIODIC JUMPING FIAS A RATIO TO WEIGHT 0.4)! "1030‘" 0 20 4.0 810 6.0 10.0 TlMEtsoc.) CONTROL POINTS 3 Aurospecrnuu E I FILE - paztnooaJnc WITH 1 91:11.; : Weight:- 186. ESTIMATE-186. lbs ‘ : .80- Val-1&1. no 11.: are l' Portodfl‘) .492 .669 .104 0 j Peak (F3) 3.222 4.636 1.396 . -40 cum (1. ) .003 n 11111941811 was ~1.026 1.00 . _ A ' t .0 2.0 4.0 6.0 ao 100 y flamumyoki Figure 4.11 CONTROLFOR Output (Periodic Jumping, 2 Hz., Vertical Ratio) PERIODIC JUMPING Fy AS A RATIO TO WEIGHT 0“!)3) "1030“ 4:0 60 TIMEuec.) AUTOSPECTHUM CONTROL POINTS FILE 8 pJ218003.FRC WITH 1 PER. Weight= 186. ESTIMATE=186. lbs Variant. we an STD '20 Period (T ) .492 .559 .104 Peak (P2) 3.222 4.536 1.395 Chain (6 ) .003 81: Peak VII .... -1.026 .00 'C”-‘DOO ‘0‘"00003 .0 O .0 20 4:0 610 6.0 10.0 firmnnoflHfl 1W1 Figure 4.10 CONTROLFOR Output (Periodic Jumping, 2 Hz., Parallel Ratio) 62 ..................................... Table 4.13 Peak Analysis (Periodic J ouncing) ......... .............. high prompt frequency. To further study this pattern, the period of each cycle of a typical periodic jumping at 2, 3, and 4 Hz. were plotted against the cycle number as shown in figures 4.12. The figure shows that when the person performs the motion at a period less than the prompt frequency, an attempt is made in the next cycle to increase the period to adjust to the prompt frequency. This assumption proved to be valid for most of the recorded cases for jumping at 2 Hz. The figure also shows that this particular person started at frequency different from the prompt. However, as the test proceeded, the person adjusted to the 63 Individual Jumping Period ( Prompt at 2, 3, and 4 Hz.) Puiod 0.5 2112.] Period (880.) 0 5 10 15 20 25 Cycle No. Figure 4.12 A Typical Periodic Jumping Response Period prompt and the variation around the prompt period decreased. On the contrary, when the prompt period decreased to 0.33 (3 Hz.) and 0.25 (4 Hz.), the two individuals were jumping consistently at lower frequencies than the prompt frequency with very few attempt to adjust to it. This observation will be used in the modeling process as explained in the next chapter. 4.4.2 Periodic J ouncing Five control points were sufficient to describe vertical force ratios for each cycle in the periodic jouncing motion One is the peak point in the cycle and the other four define the shape of the negative part. Figure 4.13 presents a typical vertical force ratio of one individual performing periodic jouncing at 2 Hz. Average peak force ratio in the vertical direction for a person jouncing in place at a 2 Hz. prompted frequency was about 0.84 times the person's weight, Table 4.15. This ratio increased to approximately 1.0 at 3 Hz. and 4 Hz. Maximum peaks also followed Vertical Force Ratios (Periodic Jouncing, 2 Hz., One Person) “W Dynanic Force Ratio b In P__.______.___——.— I...— I: .2 III: Time (Sec.) Figure 4.13 A Typical Periodic J ouncing the same pattern and increased form an average of 1.06 for 2 Hz. jouncing to 1.37 and 1.33 for 3 and 4 Hz. respectively. The average period of individuals jouncing periodically at 2 Hz. (period of 0.5 seconds) was found to be 0.501 seconds. This means that on the average, individuals were jouncing at the same prompted frequency of 2.0 Hz. As in the case of periodic jumping, individuals were unable to maintain coherency with higher frequency (3 and 4 Hz.) and their average jumping period were 0.368 (2.71 Hz.) and 0.28 (3.57 Hz.) respectively. 4.4.3 Side Swaying Side swaying forces were dominant in the horizontal direction and parallel to the individuals' seating. The recorded motion is shown in figure 4.14 where the presented force ratio is in the horizontal plane. Only peak forces were collected for the analysis. Typical force ratios in the horizontal plane are presented for the case of periodic swaying 65 1 2 3 4 5 6 7 8 9 NHHI—‘D—‘i—lh—Ip—Ip—IHH oxoooxio‘m-AuNHO Horizontal Force Ratio (Side swaying, One Individual) 1111111 ””1- :5.- Dynarnic Force Ratio o “I— "H R . I “7' vyviirygi -0.4 Figure 4.14 A Typical Swaying Side to Side 66 at 2 Hz, Table 4.16. Maximum force ratios in the horizontal direction were found not to correlate with the individuals weights. The same conclusions that were obtained from the analysis of periodic jumping and Table 4.16 Max. Peaks (Side Swaying) periodic jouncing regarding the inability of individuals to maintain coherency at higher frequencies were observed in this motion also. Table 4.17 shows typical statistics for reported average and maximum peaks in the horizontal direction. It shows the average Table 4 17 Typical Side Swaying Data force was about 0.125 times the person's , 196 0.096 0.199 weight with a range from 0.03 to 0.274. 140 0034 0-175 186 0.03 0.147 The Maximum peak ratios were in the range 150 0.053 0.149 189 0.274 0.337 of 0.112 to 0.354 times the persons' weight 125 0.197 0389 , 234 0.11 0.211 wrth an average of 0.208. 132 0.112 0.182 146 0.103 0.204 191 0.147 0.222 166 0.122 0.179 215 0.209 0.259 331 0.177 0.227 153 0.184 0.354 307 0.1 0.112 103 0.089 0.17 126 0.129 0.283 282 0.086 0.156 145 0.122 0.169 340 0.081 0.144 5.0 HUMAN LOAD MODELING 5.1 Introduction A computer program; Human Load Simulation, HLS; utilizing the understanding of the nature of each motion and factors that affect its force history was developed giving researchers an advanced tool to use in design and analysis. The forcing functions produced by the program simulate several types of human movements. The probable levels of instantaneous or peak values of human loads are represented by probability density functions. The program may be used to generate load distributions for any combination of human loading conditions, composed of those presented herein. Input to the program includes information about the type of action to be simulated, group size, and parameters necessary to define the individual loads. The program calls an array of subroutines to calculate forcing functions for different human actions and these forcing functions are combined in any form requested by the analyst. Output of the program is based on the type of action and can be summarized as follows: Transient loading (pulse loading): Statistical parameters of the control points associated with the time histories of the transient action are reported. The amplitudes of peaks and their arrival time are calculated. The program would also provide the probability density function of the peaks or the impact factors for the input action. 67 68 Steady state loads (periodic loading): The program produces force-time histories, estimates the probability densities of peaks, and computes the spectral energy distributions in the frequency domain in the form of power spectral density functions (PSD). Power spectra are determined by taking the Fourier transform of the autocorrelation functions associated with the simulated time histories, Clough and Penzien (1975). Fast Fourier Transform (FFT) method outlined by N ewland (1975), and implemented by Harichandran (1984) is adopted in the program. In order to develop the above mentioned Human Load Simulation program, HLS, several steps were required. First, model hypotheses were examined and summarized in the previous chapter. In this chapter, an attempt to describe an efficient methodology to simulate human actions is presented. The next chapter illustrates several applications to explain different ways of using the developed software. This chapter explains the random processes application to the development of the program, the generation of correlated random variables to use in estimating control points for each load type considered, and development of the Monte Carlo simulation process. The intention of the next section is to briefly explain the mathematics involved in these steps for the use in modeling human loads. 5.2 Basic Concepts 5.2.1 Random Processes Since loads due to human movements vary randomly, current methods of handling them are fundamentally deficient because they do not provide quantitative information II) 69 regarding load variabilities in a rational way. Random process theory deals directly with uncertainty by characterizing random excitations, human loads in this case, with a statistical approach. In this way, a large portion of the uncertainty in defining the forces due to human actions can be investigated and expressed in a way that can be better understood and used. The application of the concept of random processes to the modeling of the loads due to human movements seems to be very appropriate in this sense. The function describing the force due to a human action is a " random process " because it varies with time. A random process cannot be defined by a single measure or function. The set of all possible "samplefimctions" constitutes the random process, and is called the "ensemble". Figure 5.1 displays three typical samples in an ensemble of random process F(t). The figure shows that, a different sample function is obtained every time the load is recorded. The random process can then be described through its probability properties, which allow estimates or forecasts within some degree of confidence, Augusta (1984). Some of the important definitions in using random processes are described as follows: Ensemble Average: Is the average of the random process F(t), measured at time t: N F EIFUIH = ”mum; J—(ni-L .......... [5.1] i=1 where N: is the number of samples in the ensemble. Stationary Random Process: If all joint probability density functions obtained for the ensemble do not depend on time, the random process is said to be stationary. If all the average values (moments) remain constant over a change in time, the process is called 70 a "strong stationary process". If only the first and the second moments (mean and covariance) of the process remain constant with respect to time, the process is a " weak stationary process", where the first and second moments expressions are 121:17iitjii] := ‘E:[ [7( QB) ] [5:2] E[F(tl)F(t2)] = E[F(tl+T)F(t2+T)] 2 Hz. Periodic Swaying 4.14/\MAA “MIAMI/WIMP t» WVVWUVWWVWVWW.WW. inflHMMMMMM MMMN t, UUIWWWWVIVIWISWW“ ,gthWMMMMAMMAI u, “U NUQUVflVWVMWW Figure 5.1 Typical Samples in the Random Process F(t) Force Ergodic Process: The process is called "ergodic" , if in addition to being stationary, the average taken at any sample is the same as the ensemble averages. If the process is ergodic, it is stationary. 71 = Limm%,f_:ll:F(t,)dt = E(F(t))..-- [5.3] Autocorrelation: Autocorrelation is the average value of the product of the forcing function at time t and at time (t+‘c) for a random process F(t) , sampled at time t and then again at time t+t. If the process is weakly stationary, the autocorrelation depends only on t and is equal to the mean square of the process F(t): 19(1) =E[F(tl).F(tl+1:)] .......... [5.4] forergodicprocess Rf(r) = %LTF(t)F(t+r)dt, T—ooo [5.5] Power Spectral Density: It defines the frequency composition of the forcing function F(t), and is the Fourier Transformation of the autocorrelation function. Sim) = 31;]_:Rf(c)exp(-im)dt ........ [5.6] where 00 = is the circular frequency. 5.2.2 Random Number Generation Random number generation is an important component of every stochastic simulation model. An essential property of every random number generator is the ability to generate random variables that are uniformly distributed on the interval (0,1) and stochastically independent. John von Neumann (1951) and others suggested using the computer's arithmetic operations to produce sequences of numbers that, while being entirely deterministic, had the appearance of randomness. The numbers resulting are generally 72 termed pseudorandom numbers. That is, since the numbers are generated algorithmically, they are not actually random. However, for practical purposes, pseudorandom numbers are considered to behave as random numbers, i.e. , to be uniformly distributed and mutually independent provided that the random generator has a long- enough cycle length. The cycle length is the number of pseudorandom numbers generated before the same sequence of numbers starts again. Several methods may be employed to generate uniformly distributed random numbers (Motooka 1954) . The first method uses a recurrence relation to generate successive random numbers by employing any algorithm consisting of arithmetic and logical operations. The second method employs a special computer accessory, namely a physical random-number generator, that transforms the results of a physical random process (like electrical noise generators, pulse detectors of ionizing radiation events, or gas discharge tubes) into a sequence of binary digits in the computer (Golenko and Smiryagin 1960). This method increases the speed of computation, since the generated numbers are written in a fixed standard location of the memory. A third, rarely-used, method consists of inserting tables of uniformly distributed random numbers into the working memory of the computer. The fourth method, which is used in this study, is to obtain uniform random numbers by the mixing method. This method is based on utilizing specific features of the electronic computer to generate pseudorandom numbers. The method begins by assigning an initial integer number, which is placed in a certain location in computer memory. Then, the random numbers are calculated by shifting the image of the bits of the previously generated number a certain number of places. This technique of shifting 73 registers in the computer memory is used by the UNI uniform random generator subroutine to generate uniform pseudorandom numbers for this study. The advantages of pseudorandom numbers over purely random numbers is that the computer can generate them. Moreover, an important statistical advantage is the ability of such subroutines to reproduce the same sequence of pseudorandom numbers if it starts with the same initial value in the generating formula. 5.2.3 Generation of Correlated Random Numbers: From the analysis of different human load motions, it was found that correlation exists among control point arrival times and in some cases among the control point amplitudes. It is therefore necessary to preserve the same correlation structure if these variables are to be estimated as random variables to simulate the motion of an individual. The procedure explained here to generate correlated random variables was programmed by Harichandran, (1992), and is based on a methodology explained by Johnson, (1987) for the generation of multivariate normal distributions. Assume that it is desired to generate y,- random variables which maintain a correlation structure preserved in the covariance matrix cov 0’10.)- The samples of y could be generated from a sample of independent standard normal random variables x,- using the following linear transformation : [Y]=[H][X]+[p.] ............. [5.7] where H is a lower triangular transformation matrix that relates to the covariance matrix of the correlated variables y,, and 11 is a vector of the mean of each variable y,, Then it can be shown that 74 E[YYT] = E[HxxTHT] = HE[XXT]HT---- [5.8] E[YYT] = HHT ................ [5.9] since the product of E[ X X’] is an identity matrix because the random variables xi's are independent and have a unit variance. Then knowing the left hand side of eq. 5.9 which is the definition of the covariance matrix, the entries of the matrix resulting from the product H HT is found and therefore the matrix H could be obtained. Equation 5.9 represents the Cholesky factorization of the symmetric positive definite covariance matrix of Y. 5.2.4 Monte Carlo simulation Generally, a simulation model seeks to duplicate the behavior of a phenomena knowing the statistical properties of its random variables and the interactions between its components. Simulation results will reach steady state only after the numerical experiment is repeated a sufficiently large number of times. Thus, in order for the output of a model to be representative of what would be expected in the long run, simulation must produce a large enough sample (or be repeated enough times) to ensure representative results. Monte Carlo methods are those in which pr0perties of the distributions of random variables are analyzed by use of simulated random numbers. The Monte Carlo method can be defined in a broad sense as " any technique for the solution of a model using random number or pseudorandom numbers" (Hammersly and Handscamb 1964). One 75 application of the Monte Carlo method is to find the distribution or some of the parameters of the distribution of a stochastic variable (probabilistic variable). This variable is a known function of one or more other stochastic input variables that have known distributions. The method involves the determination of the distributions of the parameters, selection of a random sample of each parameter, and combination of these samples to obtain a measure of overall performance or reliability. The process of random selection and determination of response effects is repeated a large number of times, each repetition resulting in another independent estimate of the modeled function. As the sample size increases, the distribution of the sample becomes a better representation of the distribution of the population. In this study, Monte Carlo simulation is used to produce a lOOO-point sample of force-time history values for each load type combination provided by the user of HLS. 5.3 Transient Load Modeling Transient loads are modeled in HLS in a process consisting of four stages: I a) Input from the analyst provides the program with the necessary information to determine the type of transient action, number of participants, and type of output media (printer, screen, or files). I b) The program reads the parameter files for the selected motion which contains the control points statistics, phase lag parameters, and the model error (explained next). The seed array for the random number generator is read at this stage. 76 I c) HLS then calls the appropriate subroutines for the selected actions to determine the force-time history in the vertical direction, and computes the maximum positive and negative horizontal components. The model error is generated at this stage. The group force-time history is also calculated, and the maximum and minimum values of the generated forces are prepared for the output stage. I (1) Calculated individual force-time histories along with the group force time history are sent to the requested output media at this stage. The program plots force-time histories to the screen with the important statistics and saves all output files in the same directory that the program is running from. Figures 5.2 and 5.3 shows typical screen output for the case of a single jump and sitting down by one individual. In the next section, the mathematical approach used in developing load models for transient motions is explained. The single jump motion is used to demonstrate the methodology and the steps for arriving at a reliable estimate of the load produced by an individual, any number of individuals, or a group. 5.3.1 Mathematical Approach Vertical forces from transient actions are modeled in HLS as a series of lines or curves that closely fits data measured earlier and maintains the same control points statistics and correlation between different arrival times and amplitudes. Horizontal components are modeled as two impulses with magnitudes and arrival times that conforms to the observed measurements. SINGLE 77 JUMP F11 AS A RATIO TO WEIGHT VII w—v o-d)3 “030“ 58888833 28 .5 1.0 1.5 20 25 TIMEtsoc.) BIAS A RATIO TO WEIGHT 10 .oo -.10 -1 20 ’4” cl o-q)m MOIO“ 3888 .o .'6 1'0 . ‘I. TiMEtsoc.) fir 210 2.5 FZAS A RATIO TO WEIGHT 7.0:: .3 .0. I a 6.00 c 4 oo Weight- 178. lbs CONTROL POINTS Single Jump With One Person _ 1 OF 10 POINT TIES AHP. ‘ m 2 iii: ’ ' 2:23: 2.00« ‘ -' 2 1 0°. 3 .3711 1.35793 T ‘m A A, 4 .5133 -1.00699 I - VV \l’ 'V 5 .8291 -1.02749 Om-izgg \—l 6 .9285 5.16066 ’ .w I I I v 7 1.2751 “.66839 -5 1114131841651 2° 7'5 e 1. 3408 - . 00399 Romans. no (PHI) I .4337 Figure 5.2 HLS Screen Output (Single Jump, One Individual) SITTING DOWN FRAS A RATIO TO WEIGHT FyAS A RATIO T0 WEIGHT .15 .04 F I F .02 O .10 O i 1 R w. R .m C ' 0.02. E .00 ~ E ' 504 . a '1“ ' R I” 1 A A- T .610 ‘ ‘1' am ‘ ‘0 115 ' O .30 . '420 I fi I f 1 'n12 ? r .0 .5 1.0 1.5 2.0 25 .5 1.0 1.5 25 TlMEtsec.) TlMEtsec.) FZAsanArlorowmeHT CONTROL POINTS 2.60 F SITTING DOWN SIMULATION WITH 1 PER 0‘” Weight: 170. lbs '11.150 1 or 10 90121? rm: me. E 1.00 1 .0000 —.00187 2 .2785 -.51717 ’2 5° 3 .6296 .73054 T 30 4 .5083 .7723: I .50 5 .5341 1.87596 0 '00 6 .8399 -.40982 ‘1. I 1— I I l B 1.0 1.5 20 25 7 1.0113 .20364 TiMEisec.) I I Response Lag (PHI) = .4470 Figure 5 .3 HLS Screen Output (Sitting Down, One Individual) 78 Random Variables The process of transient load modeling starts by determining the weight of the person or persons performing a certain type of transient motion. Two options are available to the user of HLS ; one is to provide weights, the other to use the random number generator to estimate this variable. Currently, simulated weights are based on a model suggested by Sidney (1974). The model is a log-normal probability density function with a mean of 163 lbs, and a standard deviation of 17 lbs. The response lag 0, for any sample is predicted as a log-normal random variable with a mean of 0.27 7 seconds and a standard deviation of 0.087 seconds. Generating a log- normal random variable is based on the assumption that if a variable 0 is log-normally distributed with a mean of 11, and a standard deviation of 6,, then the natural logarithm of 0, denoted as Ln ((1)), is normally distributed. assuming )1 = 111(9)) .............. [5.10] coefficient of variation as follows: a with Coeflicient of Variation V10) = p—i’ -------- [5-11] 0 and mean and standard deviations are calculated as: 79 Estimating a log-normal random variable is a three-step process. First the transformed variable y is obtained by from Eq. [5-10] and treated as a normal random variable with a mean and standard deviation calculated from Eq. [5.12] and Eq. [5 .13]. Then a random variable is generated for the variable y , within the given mean and standard deviation, using random generations techniques discussed in Section 5 .2. Finally, the random value (1) is computed as : It is essential to note that log-normally distributed random variables cannot assume values less than zero. This assumption is suitable for the Response Lag in this case since participants always wait for the prompt to perform the intended motion. Correlated Random Variables Arrival times and amplitudes are generated as correlated random variables using the procedure in Section 5 .2. Random variables with any probability distribution function other than the normal were transformed to normal equivalence. The covariance matrix of the correlated variables were computed by multiplying the correlation coefficients by the appropriate standard deviations. The subroutine HMATRIX was used to perform linear transformation of standard normal variables to obtain the correlated variables. The mean values of the correlated variables were then added to the generated random variables. Then the new correlated random variables were transformed back to their original assumed distribution. It was found that the best-fit distributions for all random variables in the study were either normal or log-normal. Therefore, the transformation was performed using the process outlined in Eq. [5.10] through Eq. [5.14]. 80 Model Error The load model is constructed by fitting straight or curved lines to the randomly estimated control points. The determination of using straight lines or curves was based on minimizing the error or the difference between generated and measured load histories. A program MODEL.FOR, was used to fit the control points to the measured data and connect the points first by straight line and then by second degree curves. Estimates of the error associated with each trial were calculated. Figure 5 .4 is a typical screen output from the program which shows an average error calculated from 45 samples of a single jump by one individual. The error spectrum was calculated for the average error and it is shown in figure 5.5. It was found that straight lines provided less error in most cases. Figure 5 .6 shows a typical single jump where only one curve line was needed to connect control points 5 and 6. Further analysis was conducted on the error function which was defined as the difference in amplitude between measured and fitted data. Treating the error function as a random process shows that over a large number of samples, the random process approach stationarity and therefore realizations of the error can be predicted from its amplitude spectrum, Newland (1987). It is essential to note that for this particular case of single jump, a part of the process is known to be deterministic; that is the airborne time where the dynamic force ratio is always -1 . This deterministic part of the motion disturbs the assumption of stationarity. The problem could be avoided by removing the error terms from the deterministic porion of the motion. Another alternative would be to divide the error process into three M {I} III} E i... 81 ANAL‘s’SiS I3?! «'3 A 59‘ 1 :6 TC- W35“? 1... I FILE: - 33?).31344.CRL TIME 71033 (3 OP ‘4 POINT A 2 9393 - '5 . 3925‘ ‘ Q 5233 ’9 9.14"." " 6 . 9125': i 7 5. .1133 - 8 1.382;:- Rosnonu Lag (PHI) I 811' 1:117:11 1‘0 00871803.. . . ‘49:“. '80 33.111: 0.. CONTROL POINTS 5? {TE-ii 1 Weight:- 128 . ZS’E‘IPQ’I'EC 3.68 . -.1$EQ? .45313 a 05566 '.?9306 . GS???‘ .053!th . 01885 .3235 FER. £503 A”. 2331361 Ali-15.9.1465 ERROR d .00 5'; .75 7 3:- .50 . f5 . :1 .26 0‘ .w I 433‘ ..E q . Cm c‘ .'c J 1:0 126 'HMI-fizsec.) .3 in 2.5 Figure 5.4 Error Function (Single Jump, One Individual) SPECTRAL- ANKLVSIS IRI.ITf.i"$.‘-‘fil°;'l'5’.i.:M 41:H‘.1‘:i.‘.i~'11-'1.~'.I.a‘\=':t':t< A ‘3 0“! r0. 1 ‘8 ‘30. 9 1 fl .30 ' r ." a 1 Z , .20 5‘ r: ' .20 - 5 5' A 1‘. '10 j M 5 ' 1,1 v r . c H m f I i. 20.0 60.0 .0 Renews-w EH11: 210 4.0 60 30 10.0 AVEP‘ASE ~P‘ 1'."-' 9.1:“ e .18” (I -' .144 .3 .10 r .06- .02 .0 2'0 410 6.0 do 100 Freya-4113' EH 9) FRI-I: $3.11? {51644 1510843: . ‘I 2??? MEI-AK: I ii .7 MEAN: .032 1511.163 : I Sari"! {-3044 HMS: .1271" iX-Z-ZAIIK‘: $5582.10“: .132. Figure 5 .5 Error Spectrum (Sitting Down, One Individual) 82 segments; one before the airborne time, then the deterministic part followed by the landing part. The first and third segments could then be treated as two separate random processes. Removing the error term from the deterministic part was chosen in this study to simplify the modeling process. Model Error Simulation The error functions were simulated from their spectra for each generated sample and Control Polnts for a Single Jump 3 2 g. .3 1 .9 g . -1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Time (Sec.) Figure 5.6 Straight and Curved Lines of Single Jump Model then added to the simulated force ratios in the three principal axes of motion. Taking advantage of the central limit theorem and assuming that the error variable is the result of summing many statistically independent variables, the error function is assumed to approximate a Gaussian random process, Newland (1979). The mean and variance of the error stationary Gaussian process is computed for each motion considered in the study using the program MODEL.FOR. 83 To ensure consistent estimates of the error function realizations using the computed power spectrum, several smoothing band widths were considered. Figures 5 .7 and 5 . 8 are two of the error spectra where the smoothing bandwidth in the first figure was about 0.73 Hz. while it was 0.36 Hz. for the second. It can be concluded from visual inspection of the two figures that the bias introduced into the spectrum when using the higher bandwidth is not significant. The power spectrum in figure 5 .7 was used in this study to generate error functions for the case of single jump. Assuming the process to be a zero-mean, real-valued, stationary Gaussian process E(t) , and using the one-sided power spectral density function 6(0)) , the subroutine SIM U, adopted from Harichandran (1992), is used to generate samples (realizations) of E(t) using the following simulation technique, Soong (1992): l. Minimum and maximum frequencies are determined for the simulation such that the area under the one-sided power spectra could be approximated as: 5(w)=a(m),for0sws3 ......... [5.15] l The reduced power spectra is then partioned into non-overlapping intervals I: of width A0)" and amplitudes 8,, determined from the one-sided reduced spectra. I The variance 6 7,, for each interval is calculated as the area under the spectra, Figure 5 .7. l A series of independent random variables (1),, are generated to be used in the Fast Fourier Transform algorithm, FFT, to obtain a sample of size m of the estimated error function. The variable (1),, is distributed uniformly over the range (0, 21:). 84 ' The approximate sample (realization) of the error process E0) is then generated Error Spectrum (Single Jump), b=0.73 0 1 2 3 4 5 6 7 8 91011121314151617 Frequency (Hz.) Figure 5.7 Single Jump Error Function, b = 0.73 Hz. Error Spectrum (Single Jump), b = 0.36 1.00 7 0 1 2 3 4 5 6 7 8 91011121314151617 Frequency (Hz.) Figure 5.8 Single Jump Error Function, b = 0.36 Hz. 85 using the SIMU subroutine which applies a Fast Fourier Transformation on the following stochastic process, Soong (1992) to obtain the error function: Zm(t) = fiancos(wnt+()n) ......... [5.16] MI I The process is repeated for each simulation of the human load and the resulting error function is added to lines and curves connecting control points. Horizontal Force Ratios Horizontal force ratios were estimated as correlated random variables based on the measured data statistics. For each horizontal component (parallel and perpendicular to seating), the maximum and minimum force ratios were computed. Arrival times to these force ratios were determined as correlated random variables to the arrival time of the maximum force in the vertical direction. Amplitudes of the force ratios were found from the respective probability distributions, means, and standard deviations. Modeling error was established using the same process explained earlier for vertical force ratios. 5.3.2 Load Simulation The procedure explained in Section 5 . 3.1 were applied to each of the three transient motions considered in this study (single jump, sitting down, and standing up). Typical screen views of HLS output were shown earlier in figures 5 .2 and 5 .3. For each simulation, the program provides the dynamic force ratios in the three principal axes of motion along with the control points and the peak amplitudes of the sample. 86 For each motion type, HLS was used to simulate an ensemble of lOOO-samples of the motion and mean, standard deviation, and probability distributions of the control points of the ensemble were estimated. These values were then compared to the ensemble of measured data. Tables 5.1 through 5.3 present these statistics and show the absolute difference between the measured and simulated values as percentages of the measured . data. The tables demonstrate that the simulated samples converged to the measured data with differences of no more than 5 % of the measured data. 5.4 Periodic Load Modeling Periodic loads are modeled in HLS, in a process consists of five stages I a) Input to the program provides information needed to determine the type of periodic action, prompt or stimulus frequency, number of participants, and type 1.90% . Simulation 0.323 0.915 . . Tune Difference 0.20% 1.60% 1.40% 1.40% 1.30% 1.30% 1.40% 1.20% St. Dev Measured 0.086 0.041 0.106 0.117 0.192 0.187 0.262 0.286 1:12;]? Simulation 0.088 0.04 0.106 0.116 0.188 0.183 0.249 0.273 Difference 2.90% 1.30% 2.10% 5 .00% 4.40% Measured 3.155. Mean -0.062 1.134 -0.427 Aum Simulation -0.06 -O.Sl6 1.144 -0997 -0.992 3.22 -0.43 0.051 Difference 3.70% 0.50% 0.90% 0.00% 0.10% 2.00% 0.70% 0.90% St. Dev Measured 0.063 0.196 0.459 0.069 0.055 1.372 0.132 0.134 Amplit simulation 0.063 0.2 0.458 0.067 0.054 1.409 0.186 0.132 Difference 0.40% 2.00% 0.10% 2.40% 1.40% 2.70% 2.40% 1.30% 87 Measured Simulatio Difference Measured Simulatio Difference . Measured Simulatio Difference 1.40% 1.30% 0.50% 0.40% 1.40% 0.90% 1 1.10% St. Dev. Measured 0.064 0.594 0.306 0.291 0.737 0.227 0.131 Amplit. Simulatio 0.067 0.596 0.291 0.277 0.697 0.228 0.134 mffprpnm: 4 50% a 30% 4 30% o 5 40% Q 40% o of output media (printer, screen, or files). I b) The program reads parameter files which contains control points statistics, phase lag parameters, model error spectra, and the seed array for the random number generator. I c) HLS then calls the appropriate subroutines for the selected actions to determine the force-time history, and computes peak statistics for the selected motion. The model error function as a random process is generated at this stage and added to the force-time histories. I d) The one-sided power Spectra for the simulated dynamic force ratios is computed and plotted to the screen. I e) Calculated individual force-time histories along with the group force time history, if requested, are sent to the requested output media at this stage. The program is set up to plot force-time histories and spectra to the screen. Other important statistics are saved in a summary file. 88 The stochastic approach to deve10p load models for periodic motions is explained. in the next section. A first order linear Markov process is used to approximate the arrival time and amplitude peak of control points. The process as an autoregressive function is explained in detail for the periodic jumping motion. Application of the process is then extended to other motion types. 5.4.1 Stochastic Analysis Statistical analyses of measured periodic loadings revealed a consistent pattern of Table 5.3 Measured Versus Simulated Control Points For Standing up Measured Aviva! Mean Simulation 0.298 0.231 0.529 0.611 0.865 1.009 11' ms Difference 2.30% 1.00% 0.30% 0.50% 0.90% 0.60% Measured 0.17 0.12 0.181 0.195 0.37 0.379 5‘ De“ Simulation 0.173 0.118 0.175 0.192 0.377 0.379 Difference Measured Simulation 0.005 0.986 -0.672 -0.644 0.73 0.026 Difference 8.70% 0.50% 1.90% 2.20% 1.90% 1 1.10% Measured 0.016 0.33 0.35 0.39 0.268 0.088 St‘ Dev. Simulation 0.016 0.328 0.347 0.386 0.264 0.092 Difference 0.30% 0.70% 0.90% 1.10% 1.60% 5.00% Amplitudes arrival times and amplitude peaks of cycles within each recorded sample. For arrival times, It was observed that participants always attempted to adjust their motion period to match that of the prompt. This behavior was shown in Section 4.4.1 to be consistent from one individual to another and suggests some dependency between the period of the current cycle and the next one. If the probabilistic behavior at the 89 present, given all of the past behavior, depends on only the most recent past, the process is called A Markov process. Understanding the phenomena of human jumping would further suggests that the motion in a present cycle may be approximated by considering only the previous cycle. The process X is considered a first order linear Markov process, or a first order autoregressive process if it satisfies the difference equation: Xt—aXt_l = C(t) ............... [5.17] where a: is a constant and C(t) is a stationary purely random process to represent the "random disturbance" portion of the process, Priestley (1981). If the jumping period of a participant is assumed to be a first order Markov process , then Equation [3 8] could be rewritten as: T" - a TM = C ( n) .............. [5.18] where T, is the period at cycle n. Equation [5.18] could be interpreted as a linear regression equation to estimate the period T at cycle It knowing the value of the period at cycle n-1 . Tn = a Tn-l + c ( n) .............. [5.19] The process is then simplified by assuming initially conditions the period has a value of zero and Eq. [5.19] will results in : T1 = c ( 1 ) ................. [5.20] and Eq. [5 .19] is then rewritten as: 90 Tu = (n + a ( (H + a Tm2 ) ............ [5.21] Tn = (n + a (H + a2 (H + a3 (H + ,,,,,,, + art-1 T1 . . . . [5.22] and using Eq. [5 .20]: T,I = C” + a CH + a2 Cm2 + ....... + a"‘1 C1 ------- [5-231 and if E(T,_) = pg, then E(Tn) = “c(a+02+ ,,,,,,, .an-l) ......... [5.24] II EN.) = 16(1‘“ ) for a at 1 ........ [5.25] 1 — a and E(T,,) = [Lg n for a =1. If the mean of the disturbance function is zero, then the expected value of the period is zero and does not depend on the cycle number; i.e a stationary process up to the first order. If the mean of the disturbance function is not zero, then the expected value of the period will be a function of the cycle number and the process is not stationary. However, if the absolute value of the constant coefficient a is less than zero, it can be proved that for large cycle numbers n, the process T" will be " asymptotically stationary" to the first order and computed as follows, Priestley, (1981): E(Tn) = 1‘: forlargen. la|<1 ...... [5.26] A similar process is then repeated for evaluating the variance and covariance of the process T(n), and the following equations are obtained: 91 )= 5(7; T,.,) e of a for |a|<1.[5.27] cov ( T", T l-a2 11+? and it can be said that the process T, will be "asymptotically stationary" to the second order given the absolute value of a is less than one. The autocorrelation function for the process is thus: R(r)=al'l, r=0¢112131 .............. [5°30] which decays to zero exponentially as I r 1 increases. The outlined procedure from Eq. [5.17] through [5.30] is also applied to amplitude peaks which was observed to have the same pattern of the Markov property. The process of preparing and analyzing measured data consisted of the following steps: I The CONTROLFOR program was used to determine the period and peak ratio of every cycle of each sample in the ensemble of measured data. The program reports the value of the current and last period and peak ratio for each cycle in the periodic process. Table 5 .4 shows a typical sample of this process for a periodic jumping test 2 Hz for one individual. I The constant value a was estimated from each sample in the ensemble of measured data for the considered motion using correlation program CORREL. 92 Two values of the constant a are estimated from Eq. [5.28]. These are a], and a, for the arrival times T and peaks P which are the random processes of each sample. I The value of the disturbance function C(n) associated with each process was also computed and analyzed for determining its statistical properties. I The mean and the standard deviation of a, and a, were calculated over the ensemble of measured data to determine the probability distribution function that best-fit the variables. Table 5 .5 shows the statistics of the constants a, and a, for the case of periodic jumping at frequencies 2, 3, and 4 Hz. It can be seen from the values in the table that the positive correlation between each cycle and the cycle before it was stronger for the 2 Hz. jumping in comparison with the 3 and 4 Hz. jumps. This was found to be valid for both the period and the peak ratios and it also verifies the observation made earlier in chapter 4.0 about the inability of participants to adjust to the prompt at Integer frequencies higher than 2 Hz. 5.4.3 Load Modeling The periodic motions were modeled in HLS, as first order linear autoregressive processes where the period the amplitude peak of each cycle are generated from the procedure explained earlier. For each periodic motion considered in the study (Jumping, jouncing, and side swaying) the program provides the dynamic force ratios as a time series and computes the spectral density function for the simulated motion. 93 Table 5 .4 Periods and Peaks for a Typical Periodic Jumping 0.000 3.371 0.000 0.471 3.001 3.371 0.558 2.891 3.006 0.559 3.171 2.891 0.530 3.594 3.171 0.559 4.535 3.595 0.500 3.354 4.536 0.558 4.533 3.544 0.530 4.200 4.532 0.529 4.418 4.200 0.265 3.331 4.418 0.559 -1.023 3.330 0.294 3.854 4.027 0.559 3.918 3.854 0.530 4.357 3.199 0.500 2.883 4.357 0.558 3.209 2.883 0.265 -0.935 3.209 0.265 3.615 4,935 0.529 3.523 3.615 0.559 3.728 3.523 0.529 2.756 3.728 0.559 3.265 2.756 For each periodic motion, three ensembles of lOOO—samples of the motion were Simulated at the frequencies 2, 3, and 4 Hz. The mean and the standard deviation of peak force ratios and motion period were compared with measured data statistics to verify 94 Table 5.5 Constants a, and a, for the Periodic Motions the proposed model. Table 5.6 presents these statistics and shows the absolute difference between the measured and simulated values as percentages of the measured data. The table demonstrates that the simulated samples converged to the measured data with differences of no more than 1% of the mean measured data. Simulated data from all motions showed smaller standard deviations than the measured data. This is due to the fact that in measured data there were always cases where some individuals chose to end their m0tions earlier than others in the same group. Table 5 .6 Measured Versus Simulated Data (Periodic Jummng) ............................................................................................................................................... ......... 6 .0 Model Applications 6.1 Introduction This research provides an accurate model of loads due to several human actions. The model could be used to verify existing code values, provide a tool for use in reliability~ based designs, and study the effects of human loads on various special structures. Modeling Transient action would be useful in the design of elements that support mainly transient loads such as spring boards and staircases steps. Periodic loads are more likely to govern in the vibration prevention design stage especially if the loading frequency is close to the natural frequency of the structure. It has been confirmed that periodic loadings are major contributors to most of the documented cases of human discomfort due to vibration of assembly structures, Bachmann (1988). It has also been observed and verified that periodic motions would produce the maximum forces on structures when compared with transient and random motions produced by the same group of persons. It is expected that the research conclusions would provide a basic tool to be used in changing building codes to provide designs which are not only safe against catastrophic failure but serve their purpose without harmful or annoying vibrations or deflections. 6.2 Group Load Simulation A modified version of The Human Load Simulation program, HLS is used herein to generate loads for a group of participants and compare it with similar measured data. 95 96 The periodic loading module is used to illustrate the ability of the program to generate group loadings. The modified model, GROUP.FOR uses the procedure outlined in last chapter to generate periodic jumping for individuals. Force-time series from each individual is then added to a new array denoted here as the group force-time series. This process takes into account the fact that each individual would start the periodic motion at a different time (the response lag). Although the main focus in this process is to obtain the maximum force ratio of the group of individuals, the program does not attempt to align the peaks of different individuals. The assumption here is that every individual would jump independently and the overall force is the result of adding individual forces at each discrete time point in the time series. Statistics from measured data showed that the larger the group that performs the same action at the same time, the harder it is for individuals to maintain coherency and the lower the overall generated dynamic force- ratios, Table 3.7. Using the same weights from measured data, a simulation of 100- samples per group of individuals were conducted using GROUP.FOR. Output Showed that simulated mean peak forces were always lower than the measured data. Moreover, the ratio of the difference to the measured data was found to have a consistent pattern. Table 6.1 summarizes this finding. It shows the measured versus the simulated maximum peak force ratios for groups up to 40 participants. The table illustrates the possibility of a non-linear correlation between the group size and the dynamic force-ratio generated by an individual within a group. This phenomena has been addressed in earlier Stages of this project, Tuan (1981), and Mcdonald (1984). However, insufficient data hindered studying the pattern. 97 Although other potential factors may have contributed to the pattern shown in table 6.1, the " group factor " may still be the dominant part of the shown trend. It has been observed throughout the data collection process that individuals motions were always influenced by neighbors motion. Having active individuals within a group seemed to Table 6.1 Measured Versus Simulated Maximum Force Ratios energize other individuals in the same group. It is important to emphasize however, that this conclusion has not been studied statistically to justify it due to the lack of sufficient measurement. The measuring devices utilized in this study were not designed to measure individual motions of participants in a group. The output force-time histories were always the recording of total motion on the device. The program GROUP.FOR was modified to include a non-linear regression model to account for the trend explained earlier, and denoted as the " group factor ". Screen output from the modified program are shown in Figures 6.1 and 6.2. The figures show that the simulated group of 10 individuals produced a maximum peak ratio of 1.87 while the simulation for 40 individuals produced a dynamic force ratio of 0.511. Since the simulated experiment starts with a response lag that was drawn from a distribution with a small coefficient of variation, high peaks have more probability of occurring at the first portion of the force-time series. The pattern was observed in the simulated data figures 98 GROUP SIMULATION Fl AS A RATIO IO \AIEIGHT Minimum ninth/17 .. iv . v VVVIV i o-q)r moron .-" 3 8 'ZW f l r .0 20 40 $0 00 100 TIMEtsocJ SIMULATION : 1 meme... Periodic Jumping with 1 Person Weight-156. Lbs. Prompt- 2.00 Hz. ‘0‘”00'010 Jumping Frequency --> 2.002 83. Jumping Period II) .500 Sec. net. Peak Ratio --> 1.314 _- A_-_ Sinuletioa: 1 out at 10 .o 2.0 '.o 80 10.0 '.o 8 ur am To com-1mm. .. “MW” Figure 6.1 A Typical GROUP.FOR Output, 10 Individuals m GROUP SIMULATION 8838888888 ‘C"‘”¢300 FzAS A RATIO TO WEIGHT é... [LEA/1117414111411 10111.14... :. UUUVVVVVVWWWWVW‘ :E;|'V1l-Jl-J‘\:T-lc:)r\l : IE5 AUTOQWKHRUM Periodic Jumping with 40 Per(s). AVG. Wt-163.Lbs. Prompt-2.00 Hz. nub-A seeaagg Jumping Frequency II) 1.958 Hz. Jumping Period II) .511 Sec. Hex. Peek Ratio II) .626 .0 do I) 6.0 do 1110 ananufla “"V‘DOO ‘e"'flID(n Simulation: 5 out of 40 II? INTER TO 60371101.... Figure 6.2 A Typical GROUP.FOR Output, 40 Individuals '4: 99 6.1 and 6.2 reaffirm the observation Table 6.2 summarizes the modified simulated maximum peak ratios. Maximum difference between measured and simulated group loads was found to be less than 5% 6 .3 Code Verification The focus of this study was to define the force—time history produced by one individual or a group of individuals performing a predefined typical motion. Measured and simulated human loads in the research provide a limited sets of data in the much larger spectrum of activities, motions, frequencies, and participants composition. In order to transform these load values into structural design loads, other variabilities such as spatial and temporal variations must be considered. However, knowing the maximum force ranges produced from each motion simulated in the study may help in verifying current code values for assembly structures where human beings are the major source of live loads. The highest loading effects were simulated using the Human Load Simulation program, HLS developed in this study. Statistics from the collected data confirmed that 100 periodic motions produced the highest vertical load ratios. Swaying side to side in a standing position presented the maximum force ratios in the horizontal direction and parallel to the rows of seats. Maximum transverse horizontal force ratios (perpendicular to the seating) was always associated with motions that produced maximum vertical force ratios. The single jump by an individual was found to produce the highest vertical and transverse force-ratios among all actions considered in the study (Periodic and Transient). Table 6.3 is a summary of the load ratios obtained from the current research. The log- normal distribution was found to be the best-fit distribution for all the dynamic force ratios in the table. It is proposed here to consider two separate cases of loading conditions for code verifications. One case is for structures subjected mainly to the forces of one person. Typical examples would be spring boards, car seats, and steps of a stair case. The other case is assembly structures where dynamic forces are typically generated by a group of individuals. Dynamic force ratios in Table 6.3 could be utilized when designing structures from the first case. Group loading conditions will form the bases for suggesting dynamic load ratios for structures from the second case. To illustrate the use of Table 6.3, a single jump by one individual is considered. A typical individual would generate a peak vertical dynamic force ratio of about 3.2. Assuming the individual occupies an area of 3.5 square feet, UBC (1982), and has an average weight of about 163 lb, Sidney (1979), the mean dynamic force would be about 150 lb/ftl. Similarly, the mean transverse and parallel forces could be obtained from the table considering that the code assumes parallel and transverse forces to be linearly 101 distributed over a length of 28 inches. The mean dynamic forces would therefore be about 26 and 35 lb/ft for the parallel and transverse horizontal forces. The values in the vertical and the transverse directions differ significantly from the current code values of 100 lb/ft2 for the vertical force and 10 lb/ft for the transverse force. This finding is affirmed by other studies Tuan ( 1981), and McDonald (1984). It is emphasized however that the suggested values are based on statistically limited data collection. Experimental data is still needed to arrive at more reliable estimates of load ratios. Periodic motions would be utilized for load verification of the second type of structures; assembly structures where dynamic loads are typically produced by a group of individuals. It has been found that coherency among individuals performing transient actions is generally low and therefore it significantly reduces the overall group force- ratio. Periodic jumping at about 2 Hz. was found to produce the largest vertical and transverse force ratios. The inability of individuals to maintain perfect coherency tends to reduce the group force ratios when more than one person attempts the action. It has been observed however, that a group effect exists in periodic motions where individuals may influence the motions of others in the same group. Taking both observations into account and reviewing the measured and simulated data, it is proposed to use the "Dynamic group factor" introduced in chapter 3 .0 to modify the mean dynamic force ratios in Table 6.3. The Dynamic Group Factor, denoted as DF, has the following form: DF = 0,97" 4 Fr ................. [6.1] 102 1",: Mean dynamic force ratio of one individual n : Crowd size Equation [6.1] should be used to reduce the mean dynamic force ratio for periodic actions in Table 6. 3. The corrected force ratios could then be employed to estimate an upper bound to the dynamic force by a group of individuals performing any of the actions simulated in the study. When periodic jumping was used to verify code values, it was found that the code values would provide acceptable maximum ranges for both the vertical and the parallel components of the periodic motions. However, the transverse component was found to be significantly higher than the code value. It is therefore suggested to increase the load value of 10 lb/ ft in the transverse direction to 45 lb/ ft. Other literature showed the finding to be valid. Tuan (1983), suggested a value of 43 lb/ ft for the same force component. 6.4 Reliability-Based Design Generally, the reliability of any system is defined to be the probability of performing adequately for the period of time intended under the operating conditions encountered. The reliability of a structure is also defined as the probability that the structure will not exceed any of its limit states during a specified time period. In the case of human, one of the limit states loadings on assembly structures could be considered as the probability of not exceeding a certain vibration level. A certain value of acceleration may be another form of a limit state. Defining the statistical parameters of the dynamic human loads is therefore an important and essential ingredient in finding the overall reliability of the 103 Table 6.3 Dynamic Force Ratios for One Individual structure under analysis. A currently used measure to describe structural reliability is the reliability index, B. It represents the distance, in Standard deviations, from the set of variable values that is most likely to occur to the set of variable values that is most likely to cause failure. The reliability index 13 is directly pr0portional to the difference between the mean values of the load L and resistance R. where 6L, O'R are the standard deviations of load and resistance respectively. It should be noted that this formulation for the reliability index 104 is valid only if both of the load and resistance are normally distributed. Moreover, the same principles apply for other distributions. R -L 13 - Z "'2 .................. [6.2] o, + 01. Although the reliability index does not directly indicate a value for the reliability of a structure, it can be used to indicate a lower bound on the reliability and also to serve as a relative measure for the level of safety of different structures. Beyond that, we can only state that, for a given structure, larger values of the reliability index represent greater reliabilities than smaller values. The reliability index usually varies from 2.0 to 8.0. Analysis of 3 values for specifications governing the design of buildings in the United States indicates that B typically ranges from 3.0 to 4.0. These values are for members designed to resist gravity loads in situations where failures are ductile, Ellingwood (1980). It is however important to emphasize at this point that, whereas the reliability index (B) is not a direct measure of the reliability, it is a useful comparative measure and can serve to evaluate the relative safety of various structures. The advantage in simulating dynamic loads produced by human actions may greatly benefit the reliability design concepts. Knowing the load statistics enables a designer to simulate over all probable values for loads and resistance and then determine the acceptable level of safety. The load ratios developed in this study for several human movements could be used in design cases where the knowledge of load variabilities would 105 enhance the ability to assume the appropriate safety measures and reliability for the intended structure. 6.3 Other Potential Applications It is anticipated that the results of this research will be very useful for a broad range of building industry. following are some potential examples: Floor Vibrations Modern structures which have slender and elegant elements can be very sensitive to vibrations. Increasing numbers of excessive floor vibration incidents caused by human movements have been reported recently, Buchmann (1992), and Sterareh (1992). The coincidence of the natural frequencies of the floor system with human movement frequency is the major cause of this problem. The correspondence of the structure' 8 natural frequency with the frequency of a higher harmonic of the time history of the dynamic excitation has also been realized to lead to resonance. To reduce the effect of annoying vibrations on the serviceability of structures, several solutions have been suggested. Among them, are the modifications of the floor stiffeners or mass, increasing the system damping, and/ or limiting the vibration amplitudes. The accurate representation of the dynamic forces applied on the structure is essential to develop a reliable solution to the excessive vibration problem. Manufacturers of floor systems, timber, pre-stressed concrete, open-web steel joists, and other combinations of geometry and materials may be willing to conduct a computer 106 simulation to anticipate potential problems rather than doing a full scale test or risk an adverse design for prototype installation. Safety Studies Insurance companies concerned with adverse or potentially harmful loading due to human actions or occupancy, such as occurred at the disastrous Kansas City Hyatt sky bridges, could be interested in a tool which has the capability of simulating a situation in a computer prior to construction, or in another scenario, analyzing an occurrence so as to remedy a problem. Education The Human Load Simulation program, HLS, could be used as an educational tool for graduate courses dealing with reliability design, random processes, theory of random vibrations, and dynamic of structures. The ability of the program to simulate several loading conditions (transient and Periodic) and to produce pseudorandom processes could help in illustrating many of the concepts and applications of these courses. An instructor could utilize the random number generation features in HLS to produce random variables, correlated variables, or Markov random processes for any application since the produced variables would be dimensionless. 7 .0 CONCLUSIONS 7 .1 Summary In addition to many other possible applications, human loads are the primary concern in the design and analysis of assembly structures. Human loads vary from event to event in an apparent pattern and they can be shown to be predictable statistically. The focus of this research is to quantify the uncertainty associated with the loads due to human movements. Providing accurate models of probabilistic functions to represent forces due to several predefined human movements was the main objective of the study. A detailed review of current state of knowledge in the field of modeling human loads is presented. The study also provides an overview of current practice in the design process and how the present U.S. codes characterize human loads on assembly structures. Data taken from tests performed on a force platform and a floor system were obtained from three sets of experiments; two designed to measure the force due to the action of an individual and the third designed to measure the force due to the action of a group. Two different test devices were utilized. A 4 ft by 8 ft force platform was used for measuring the vertical and horizontal components of loads produced by individuals and small groups up to five participants. The other device was a 12 ft by 15 ft floor system designed to measure total vertical loads generated by a group of people of up to forty participants. The study includes six different in-place human movements grouped into two categories; transient motions which includes single jumps, standing up suddenly 108 from a sitting position, and dropping into a seat from a standing position; and periodic motions which includes periodic jumping, jouncing, and swaying side to Side. For each transient action in the study, a method based on passing straight or curved lines through control points that define the force-time history is used to simulate the forces due to human actions. Best-fit distributions for the control points were determined. An error function was developed and incorporated into the model to better represent the modeled actions. Periodic loadings were analyzed as a series of impulses where the forcing function is divided into cycles each of which is a full period of motion considered. Statistics of each impulse in the force-time history were computed and analyzed for correlation and peak distribution. The motions were modeled as first order autoregressive processes. An iterative model-building methodology whereby the stochastic models are related to actual time series data was conducted. Model verification was applied to the different motions considered in the study. A model to simulate group loading is provided. The developed forcing functions are compared with experimental data and feed back from the process is used to refine the forcing functions. A computer program, Human Load Simulation, HLS; was developed giving researchers an advanced tool to use in design and analysis. The forcing functions developed by the program simulate several types of human movements. The program may be used to simplify generating load distributions for any combination of human loading conditions, similar to those presented in this research. Input to the program includes information about the type of action to be simulated, group size, and parameters 109 necessary to define the individual loads. The program calls an array of subroutines to calculate forcing functions for different human actions and these forcing functions are combined in any form requested by the analyst. Output of the program is based on the type of action and includes the following: I Statistical parameters of the control points. I Amplitudes of peaks and their arrival times. I Spectral energy distributions for periodic motions. I Force-time history for individuals. I Group force histories (for group loading). Several potential applications of the program developed from the study are presented. Code values were examined and recommendations regarding the transverse horizontal force is outlined. Other potential applications of the research outcome are also discussed. 7.2 Conclusions This study indicates several conclusions with respect to the modeling of human load movements: I Measured data Showed that human movements can be reproduced using probabilistic procedures. I Periodic motions produced the highest vertical load ratios. I Swaying side to side in a standing position presented the maximum force ratios in the horizontal direction (parallel to the seatings). 7.3 110 Maximum transverse horizontal force ratios (perpendicular to the seatings) were always associated with motions that produced maximum vertical force ratios. The single jump was found to produce the highest vertical and transverse force ratios among the actions considered in the study (periodic and transient). Log-normal distribution was found to be the best-fit distribution for all peak force ratios. The developed program HLS, will simplify generating load distributions for any combination of human loading conditions, similar to those presented in this research. Code values were examined and found to be acceptable except the transverse force component where a value of 45 lb/ ft is recommended to replace the current 10 lb/ ft value. Future Work The study leads to the following areas of needed further work: I The human load program developed in the study needs to be expanded to include other types of human loadings. I More statistical work is needed to establish the model error for human loads provided by HLS. I A study is needed to understand the effect of area of occupancy, floor type, and coherency among participants on the overall force ratios of one individual in a 111 group. A state of the art force platform would enable measuring individual forcing functions for individuals in a group. Aerobic motions are interesting class of human motions which deserves to be studied and modeled. Higher jumping frequencies and higher degree of coherency are expected in this class of motions. Studies on the effect of the prompt frequencies on the output forcing functions is still needed. Available data are very limited in this perspective. Taking advantage of the developed software may raise other potential possibilities of using CAD interface systems to show the performed motions and generated forces in real-time. Additional work is needed to translate the results of this research into formats suitable for use in design specifications. APPENDIX A DATA PROCESSING PROGRAMS lldinthZhR $STORAGE:2 $LARGE* * * * * * * * * * * * * * * * * * * * * i * * * * * * i * * * * * * g. ‘- 4 It I- 4 :1- :1- :1- 1(- * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 44 44 44 44 4* 44 44 44 APRIL 15, 1992 MAY 05, 1992 June 15, 1992 June 20, 1992 BASSEM K. KHAFAGI AUG 16, 1992 Aug 24, 1992 Aug 30, 1992 Oct 5, 1992 MICHIGAN STATE UNIVERSITY EAST LANSING, MI 48823 U.S.A * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * i * ***************************************************************** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ *4-44'4*=$*fl-*N-**I+*X-**ik*flr** *fi-**'**ik*fl'**it*I-#*'441+41'44 0000004444444444444444444444» PROGRAM DATA.FOR DIMENSION XVOLT(512,11), FX(512), FY(512), FZ(512) DIMENSION BVOLT(1,11), VOLT(512,11), PROMPT(512) DIMENSION TIME(512), PR(512), PRP(200) DIMENSION XCALIB(11),YCALIB(11),ZCALIB(11) DIMENSION XCALIBA(11),YCALIBA(11),ZCALIBA(11) DIMENSION XCALIBB(11),YCALIBB(11),ZCALIBB(11) DIMENSION FXP(512) FYP(512), FZP(512) DIMENSION Fme(2001, Fme(200), Fme(200),twgt(200) INTEGER WEIGHT(5),NAA(20) LOGICAL*1 XABS CHARACTER REPORT*12,INPUTS*40,JOB*40,dataset*1,TT*5,EXT*12,DP*1 CHARACTER FNAME1*3,FNAME2*3,FNAME3*3,FNAME4*3,FNAME*12,0UTPUTS*40 CHARACTER Filefla *2,fileflagn*2 CHARACTER*8 NAME( 0,200 COMMON /CMP/ M,N,FS,IWIN,L,NFFT,DF COMMON /FILES/ INPUTS,OUTPUTS,LINP,LOUT FLAG FOR THE OPEN SCREEN OOO DP='y' OFLAG=1 Read INPUT Data from PARAMETER FILE "PARAM.DAT" COO OVJA '11 E El 112 113 DETERMINE THE LOCATION OF INPUT AND OUTPUT FILES 1.0 INPUTS OOOOOO LINP= LEN TRIM(INPUTS) IF(LINP. NE. 0) THEN INPUTS= —INPUTS(1: LINP)//' \' LINP=LINP+1 END IF 2.0 OUTPUTS LOUT=LEN TRIM(OUTPUTS) IF(LOUT. NE. 0) THEN OUTPUTS=OUTPUTS(1: LOUT)//' \' LOUT= LOUT+1 END IF START FILES' PROCESSING READ (1,*) NG DO 5 LL=1,NG READ (1,*) NA NAA(LL)=NA DO 20 I= 1, NA READ (1,16) FNAMEl, FNAMEZ, FNAME3, DATASET, FNAME4 10 FORMAT (A2, A1, A1,A1,.A3) 20 WRITE (NAME(LL, I), 30) ENAME1,ENAME2,ENAME3,DATASET, PNAME4 30 FORMAT (A2, A1 AL A1, A3) IF (LL. EQ. 1 THEN WRITE (EXT,9)JOB,fname2 GOO 000 9 FORMAT (A2,'-',a1,'-SUM.OUT') OPEN (7,FILE=OUTPUTS(1:LOUT)//EXT) END IF 5 CONTINUE C C C Read spectral estimation parameters C READ (1,*) M READ (1,*) IWIN READ (l,*) L READ (1,*) NFFT READ (l, *) OF CLOSE (l) C C c +++++++++++++++++++++++++++++++++++++++ C +++++++++++++ 2. READING CALIBRATION EACTORS +++++++++++++ C +++++++++++++++++++++++++++++++++++++++ C C C READ THE CALIBRATION FACTORS FILES IN THE THREE DIRECIONS X,Y,Z C USE THESE FACTORS TO CALCULATE THE FORCES IN EACH DIRECTION BY C MULTIPLYING THE VOLT(I,J) * CALIB-X FOR THE X FORCES AND THE SAM C FOR THE OTHER TWO DIRECTIONS Y,Z C C READ CALIBRATION FACTORS FOR THE X, Y, Z FOR DATA SET A,B C C C OPEN THE CALIBRATION FILES IN THE X, Y, Z DIRECTIONS. C C DO 35 NN=1,2 IF (NN.EQ.1) THEN OPEN (3,FILE='X-CALIB.DSl') OPEN (4,FILE='Y-CALIB.DSl') OPEN (5,FILE='Z-CALIB.D81') DO 40 J=1,11 READ (3,50) XCALIBA(J) READ (4,50) YCALIBA(J) READ (5,50) ZCALIBA(J) 40 CONTINUE CLOSE (3) CLOSE (4) 45 35 50 000 000 60 70 000 75 85 000 80 0000 90 100 110 00000 114 CLOSE (5) ELSEIF (NN.EQ.2) THEN OPEN (3,FILE='X-CALIB.DSZ') OPEN (4,FILE='Y-CALIB.DSZ') OPEN (5,FILE='Z-CALIB.DSZ') DO 45 J=1,11 READ 3,50) XCALIBB(J) READ 4,50) YCALIBB(J) READ (5,50) ZCALIBB(J) CONTINUE CLOSE (3) CLOSE (4) CLOSE (5) END IF CONTINUE FORMAT (11(E18.3)) START THE FILES-PROCESSING LOOP DO 480 LL=1,NG NA=NAA (LL ) DO 380 K=l,NA E.G. THE LOOP CAN BE CHANGED TO K=1,5 IF WE WANT TO PROCESS 5 FILE WRITE (FNAME, 60) NAME(LL, K) FORMAT (A8, W) WRITE (REPORT, 70) FNAME FORMAT (A8,'.FRC') READ (FNAME,10) FNAMEl,FNAME2,FNAME3,DATASET,FNAME4 determine the calibration files to use and the sampling rate for the test IF é;DATASET.EQ.'A}).OR.(DATASET.EQ.'a')) THEN FS= . DO 75 J=l,11 XCALIB(J)=XCALIBA(J) YCALIB(J)=YCALIBA(J) ZCALIB(J)=ZCALIBA(J) CONTINUE ELSEIF ((DATASET.EQ.'B').OR.(DATASET.EQ.'b')) THEN FS= . DO 85 J=1,11 XCALIB(J)=XCALIBB(J) YCALIB(J)=YCALIBB(J) ZCALIB(J)=ZCALIBB(J) CONTINUE ENDIF BEGIN BY OPENING FILES TO READ DATA OPEN (2, FILE= INPUTS(1: LINP)//FNAME) OPEN (6, FILE= OUTPUTS(1: LOUT)//REPORT) READ (2,80) (WEIGHT(I), I= 1, 5) FORMAT (/, 5(I4)//) Now add the weights of participants to get TWEIGHT TWEIGHT= O. 0 DO 90 I=1,5 TWEIGHT= TWEIGHT+WEIGHT(I) CONTINUE IF((Fname2.eq. 't'). or. (FnameZ. eq. 'T' ). or. (Fname2. eq. 'R' ). OR. 1(Fname2.eq.'r')) N= 256 D0 100 I=1,N READ (2,110) (XVOLT(I,J),J=l,11),PROMPT(I) CONTINUE FORMAT (11(F8.0/),F6.0) +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 3. NORMALIZING VOLT. READINGDS +++++++++++++ +++++++++++++++++++++++++++++++++++++++ 000000 120 00000 000 144 134 000 140 130 000000 150 160 0000 000 115 SELECT THE FIRST (OR THIRD) VOLTAGE READING AS THE BASE TO NORMRLIZE THE VOLTAGE READING AGAINST THE PART. WEIGHTS AND SAVE IT IN MATRIX "BVOLT" DO 120 I=1,11 IF ((DATASET.EQ.PA').OR.(DATASET.EQ.'a')) THEN BVOLT(1,I)=XVOLT(1,I) ELSE BVOLT(1,I)=XVOLT(3,I) END IF CONTINUE DETERMINE THE TIME STEP USED IN COLLECTING DATA TIMESTEP=1.0/FS NORMALIZE DATA BY SUBTRACTING EACH RAW OF DATA FROM THE FIRST VLOTAGE READING WHICH REPRESENT THE STATIC WEIGHT OF THE FLOOR AND PARTICIPANTS. DO 134 I=1,N DO 144 J=1,ll VOLT(I J)=0.0 TIME(I)=0.0 PR(I)=0.0 CONTINUE CONTINUE DO 130 I=1,N DO 140 J=1,11 VOLT(I J)=XVOLT(I,J)-BVOLT(1,J) TIME(I)=TIME(I—l)+TIMESTEP PR(I)=ABS(PROMPT(I)-PROMPT(1)) IF (I.EQ.1) THEN TIME(I)=0.0 END IF CONTINUE CONTINUE +++++++++++++++++++++++++++++++++++++++ +++++++++++ 4. FINDING FINAL FORCES X, Y, z , P +++++++++++ +++++++++++++++++++++++++++++++++++++++ DO 150 I=1,N FX(I)=0.0 FY(I) FZ(I) FXP(I FYP(I FZP(I CONTINUE DO 170 I=1,N DO 160 J=l,11 FX (I)=FX (I) +VOLT (I, J) *XCALIB(J) FY(I)=FY(I)+VOLT(I,J)*YCALIB(J) FZ(I)=FZ(I)+VOLT(I,J)*ZCALIB(J) CONTINUE SWITCH FX WITH FY FOR DATA SET A (TO MATCH DATA SET B COORD. SYSTEM) aLSO TAKE THE WEIGHT OFF FROM THE FZ(I) FOR DATA.SET A .0 .0 II II ll CO 000 .0 .0 .0 IF ((DATASET.EQ.'A}).OR.(DATASET.EQ.'a')) THEN FXSWgFX(I) FYSW=FY(I) FX(I)=FYSW FY(I)=FXSW FZ(I)=FZ(I)-TWEIGHT ENDIF FIND FORCE RATIO 170 000000 000000000000000 000 180 0000 116 FXP(I)=FX(I)/TWEIGHT FYP(I)=FY(I)/TWEIGHT FZP(I)= Z(I)/TWEIGHT CONTINUE NOW FIND THE MAXIMUM VALUE OF THE FORCE FX,FY,FZ, AND USE EACH AS A.VALUE FOR THE PROMPT SPIKE. CALL MAX (FX,FXMAX,TX,TIMESTEP) CALL MAX (FY,FYMAX,TY,TIMESTEP) CALL MAX (FZ,FZMAX,TZ,TIMESTEP) CALL MAX (FXP,FXPMAX,TX,TIMESTEP) CALL MAX (FYP,FYPMAX,TY,TIMESTEP) CALL MAX (FZP,FZPMAX,TZ,TIMESTEP) CALL MIN (FX,FXMIN,TXM,TIMESTEP) CALL MIN (FY,FYMIN,TYM,TIMESTEP) CALL MIN (FZ,FZMIN,TZM,TIMESTEP) CALL MIN (FXP,FXPMIN,TXM,TIMESTEP) CALL MIN (FYP,FYPMIN,TYM,TIMESTEP) CALL MIN (FZP,FZPMIN,TZM,TIMESTEP) FOR THE CASE OF SINGLE JUMP OR PERIODIEC JUMPING WHERE THE PERSON LEAVE THE FLOOR, IT MAY BE POSIBLE TO ESTIMATE THE STATIC WEIGHT OF THE PERSON(S). THIS IS USEFUL IF THE WEIGHT IS UNKNOW OR AS AS MEASURE OF DATA AQUISITION ACCURACY. IT HAS BEEN FOUND THAT A PERSON LEAVES THE FLOOR WHEN JUMPING FOR ABOUT 0.6 SECONDS (ABOUT 20 READINGS) AND COMES FOLLOWED DIRECTLY BY THE FZMAX (THE HIT). EWEIGHT=0.0 IF ((FNAME1.EQ.'SJ').OR.(FNAMEl.EQ.'sj').OR.(FNAMEl.EQ.'PJ').OR. 1(FNAME1.EQ. 'pj ' ) ) THEN EWEIGHT=-fzm1n END IF REPORT THE LOCATION OF THE PROMPT (IF ANY) FOR TRANSIENT ACTION IF ((ENAME2.EQ.'t').OR.(ENAME2.EQ.'T')) THEN CALL MAX (PR,PRMAX,TP,TIMESTEP) DO 180 I=1,N PR(I)=0.0 CONTINUE NP=INT(TP*FS)+1 PRMAX=100. PR(NP)=PRMAX ELSE FIND PEAKS FOR THE PROMPT FOR PERIODIC ACTION CALL MAX (PR,PRMAX,TP,TIMESTEP) CUTOFF=PRMAX/4.0 NPK=O DO 190 I=1,N 190 000 200 210 000000000 220 230 240 00000 250 00000 260 270 000 280 000 290 117 IF (PR(I).LE.CUTOFF) THEN PR(I)=0.0 ELSE NPK=NPK+1 PRP(NPK)=I PRMAX=100. PR(I)=PRMAX END IF CONTINUE THIS LOOP IS TO ENSURE ONLY ONE READING PER PROMPT DO 210 Nz=1,NPK I=PRP(Nz) DO 200 NN=1,4 PR(I+NN)=0.0 CONTINUE CONTINUE END IF +++++++++++++++++++++++++++++++++++++++ +++++++++++++ S. START THE FINAL OUTPUT +++++++++++++ +++++++++++++++++++++++++++++++++++++++ PRINT TIME VS FORCES IN X,Y,Z DIR. PLUS THE PROMPT WRITE (6,220) , ' FORMAT (///,10X,'TIME ',3X,'X-FORCE Ratlo ',2X,'Y-FORCE Ratlo', 12X,'Z-FORCE Ratio',2X,'PROMPT',/) DO 230 I=1,N WRITE (6,240) TIME(I),FXP(I),FYP(I),FZP(I),PR(I) CONTINUE FORMAT (7X,F8.3,3(F14.6),3X,F10.0) \ PRINT INFORMATION ABOUT THE FILE NAME TO THE SCREEN TO ENSURE THAT CORRECT FILE IS BEING HANDLED WRITE (6,250) FNAME,FNAME3 FORMAT (//,20x,' INPUT FILE INFORMATION',///,3x,'FILE NAME = ', 1A12,10x,'PARTICIPANTS = ',A3,2x,'PERSONs.'/) REPORT AN EXPLAINATION OF THE INFORMATION STORED IN THE FILE THIS INFORMATION WILL BE PRINTED TO THE SCREEN AND SAVED IN THE REPORT FILE (*.#RP) IF ((FNAME2.EQ.'t').OR.(FNAME2.EQ.'T')) THEN WRITE (6,260) FNAME1,FNAME4 ELSE WRITE (6,270) FNAME1,FNAME4 END IF FORMAT (3x,'TRANSIENT ACTION FORMAT (3X,'PERIODIC ACTION THEN, PRINT OUTPUT SUMMARY WRITE (6,280) (WEIGHT(I),I=1,5) TWEIGHT,EWEIGHT FORMAT }//,1X,'WEIGHT(S) z',5(Ié),3x,'Tota1 :',F7.0,5X,‘ Estimate ) ',A3,5x,'TEST NO. = ',A3/) ',A3,5x,'TEST NO. = ',A3/) 1:',F7.2 FINALLY, WRITE A SUMMARY ABOUT THE FORCES. WRITE (6,290) REPORT FORMAT (20x,' FORCES INFORMATION ',///,13x,'FINAL REPORT FILE N 1AME = ',A12///) WRITE (6,300) FXMAX,FXPMAX,TX WRITE (6,310) FYMAX,FYPMAX,TY WRITE (6,320) FZMAX,FZPMAX,TZ IF ((FNAME2.EQ.'t').OR.(FNAME2.EQ.'T')) THEN WRITE (6,330) PRMAX,TP ELSE WRITE (6,340) PRMAX END IF 000000000000000 0.30000 0000000000 000000000 300 310 320 330 340 50 118 FORMAT (2x,'FXMAx =',F12.4,' LBS RATIO TO WEIGHT=‘,F8.4,8X,'AT T 1IME = ',F5.3,' SEC.',/) FORMAT (2X,'FYMAX =',F12.4,' LBS RATIO TO WEIGHT=',F8.4,8X,'AT T lIME = ',F5.3,' SEC.',/) FORMAT (2X,'FZMAX =',F12.4,' LBS RATIO To WEIGHT=‘,F8.4,8X,'AT T lIME = ',F5.3,' SEC.',/) FORMAT (2X,'PRMAX =',F12.2,' ',26x,'AT TIME = ',F5.3,' SEC.',/) FORMAT (2X,'PRMAX =',F12.2,' ',26x,'AT TIME = PERIODIC',/) +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 6. PREPARING PLOTS TO THE SCREEN +++++++++++++ +++++++++++++++++++++++++++++++++++++++ CALL THE OPENING MENU FOR FINAL ON—SCREEN OUTPUT(UNDER DESIGN) THE OPEN FLAG‘WILL OPEN THE MENU ONLY THE FIRST TIME THE PROGRAM THE COLOR FLAG WILL BE SET TO (0) OR (1) DEPENDING ON THE ITERATION NO. THIS FLAG'WILL INTERCHANGE THE BACKGROUND COLOR. IF (OFLAG.GT.0.0) THEN IF((DP.EQ.'Y').OR.(DP.EQ.'y')) THEN CALL OPEN (OFLAG,BW) READ (*,*) END IF ELSE GO TO 350 END IF CALL PLOTXY SUBROUTINE TO PLOT THE TIME VERSUS THE FORCE IN THE X, Y, z DIRECTIONS. CONTINUE IF((DP.EQ.'Y').OR.(DP.EQ.'y')) THEN CALL PLOTXY (TIME,FXP,FYP,FZP,PR,FXPMAX,FYPMAX,FZPMAX,NP,TX,TY,TZ, 1WEIGHT,TWEIGHT,CFALG,K,NA,FXPMIN,FYPMIN,FzPMIN,TXM,TYM,N,FNAME, 2EWEIGHT,BW) ELSE print *,'GROUP: ',ll,'FILE: ',k END IF CLOSE (2) CLOSE (6) ++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++ 7. SPECTRAL ANALYSIS ESTIMATION OF FORCE-SERIES ++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++ spectral analysis will be called out only if the action is periodic. Also if more than on action is being processed, the average correlation and average spectrs Will be disabled IF((Fname1.eq.'ss').or.(Fname1.eq.'SS').or.(Fnamel.eq.'jn').or. l(Fnamel.eg.'JN').or.(Fname1.eq.'PJ').or.(Fnamel.eq.'pj')) THEN if(k.eq. ) then noav%=0 file lag=fname1 endif fileflagn=fname1 if(filef1agn.ne.fileflag)noavg=l CALL SPECTRA (LL,K,NA,NAME,OF,RMS,ZMEAN,BW,noavg) END IF ++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++ 8. Reporting max forces to file "SUMMARY.OUT" ++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++ Find absolute maximum forces and report them.to the matix fxmx(i) fxmx(k)=fx max . if(fxpmax. t.abs(fxpmin)) fxmx(k)=abs(fxpmin) fymx(k)=fypmax 000 00000 000000 0000) 000 000 000 360 370 380 371 372 480 119 if(fyim ax. 1t. abs(fypmin)) fymx(k)=abs(fypmin) fzmx ) =fzpmax twgt(k)=tweight IF (K.EQ.1) THEN WRITE (7,360)fname3,FNAMEZ ENDIF WRITE (7, 370) K, TWGT(K), FXMX(K), FYMM(K), FZMM(K) FORMAT(//, 5x, SUMMARY OF TESTS OF' ,As, PERS. AT ',al,‘ Hz.',//, llX, 'Test 'Weight Fxmax Fymax Fzmax') FORMAT (2X,IZ, 5X, F8. 2, 3(4X,F6.3)) DONE ..... print *,ll,k CONTINUE CALCULATE AND REPORT IMPORTANT STATISTICS FOR THE OBSERVED ACTION This is done only if the same action is considered if(noavg.eq.0)then XABS=.TRUE. CALL MYSTAT (TWGT,NA,WMIN,WMAX,NWIN,NWAX,WMEAN,WVAR,WSTD,WVX, 1W5KE,WKUR/WSE,XABS) CALL MYSTAT (FXMH,NA,XMIN,XMAx,NXIN,NXAX,XMEAN,XVAR,XSTD,xvx, 1XSKE,XKUR,XSE,XABS) CALL MYSTAT (FYMX,NA,YMIN,YMAX,NYIN,NYAX,YMEAN,YVAR,YSTD,YVX, 1YSKE,YKUR,YSE,XABS) CALL MYSTAT (FZMx,NA,2MIN,ZMAx,NzIN,NZAx,ZMEAN,ZVAR,ZSTD,zvx, lZSKE,ZKUR,ZSE,XABS) BLANK LINE T0 SEPARATE DATA WRITE (7,373) FORMAT (1X) REPORT TO THE OUTPUT FILE TT='MEAN' WRITE (7, 371) TT, WMEAN, XMEAN, YMEAN, ZMEAN FORMAT (1X,.A5, 1X, F10. 2, 3(2X, F8. 3)) FORMAT (1X, A5, 1X, IlO, 3(2X, I8), 3X, Ill) TT='MIN' WRITE (7, 371) TT,WMIN,XMIN,YMIN,ZMIN TT=' MAX' ¥¥ITE (7, 371) TT,WMAX,XMAX,YMAX,ZMAX WRITE (7, 371) TT,WVAR,XVAR,YVAR,ZVAR TT='STD' WRITE (7,371) TT,WSTD,XSTD,YSTD,ZSTD TT='C.O.V' WRITE (7,371) TT,WVX,XVX,YVX,ZVX TT='NMIN' WRITE (7, 372) TT,NWIN,NXIN,NYIN,NZIN TT=' NMAX' WRITE (7,372) TT,NWAX,NXAX,NYAX,NZAX TT='SKEWN' WRITE (7,371) TT,WSKE,XSKE,YSKE,ZSKE TT='KURTO' WRITE (7,371) TT,WKUR,XKUR,YKUR,ZKUR TT='ST- -ER' WRITE (7, 371) TT,WSE,XSE,YSE,ZSE end CONTINUE CLOSE(7) ALL DONE CALL ENDGRAPHICS() 120 STOP END MAN PROGRAM DATArP #################################f##################################### * * k * * * * * * * * * * * * * * * * * * * i * * * * * * * * * * * x» ** *- SUBROUTINE OPEN * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * APRIL 15,1992 OPENING LOGO...... * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ***************************************************************** **-#*-**’*W-** *fi-*i-*1~*i-** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ 0000000***w****»****00 SUBROUTINE OPEN (OFLAG,BW) LOGICAL*1 INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, 1PLOTMETHOD,INVERTX,INVERTY,TITLEONLY,SCRONLY,Y2PLOT,y3plot, 2PLOT2METHOD INTEGER*2 BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 1TITLECOLOR,NUMBER,PLOTZCOLOR,FONTHEIGHT,FONTWIDTH,PLOT3COLOR CHARACTER TITLE*280,XAXISTITLE*80,YAXISTITLE*80 C INCLUDE 'PLOT.DIF' C 7 INITIALISE THINGS FOR THE FIRST PLOT INITGRAPHICS=.TRUE. INITBORDER=.TRUE. INITWINDOw=.TRUE. INITAXES=.FALSE. INITFONTS=.TRUE. PLOTMETHOD=.TRUE. Y2PLOT=.FALSE. Y3PLOT=.FALSE. SCRONLY=.TRUE. FONTHEIGHT=28 FONTWIDTH=16 C SET UP THE FANCY COLORS BORDERCOLOR=14 WINDOWCOLOR=15 IF (BW.EQ.1) THEN BORDERCOLOR=0 WINDOWCOLOR=15 END IF C SET UP THE PLOT LIMITS XMIN=0.0 XMAX=1.0 YMIN=0.0 YMAX=1.0 C POSITION THIS PLOT AT THE TOP LEFT OF THE SCREEN XBOTTOMLEFT=0.0 YBOTTOMLEFT=0.0 XTOPRIGHT=1.0 YTOPRIGHT=1.0 TITLE=' ** D A. T .A . F O R **' CALL PLOT (INITGRAPHICS, INITBORDER, INITWINDOW, INITAXES, INITFONTS, 1XDATA, YDATA, Y2DATA,NUMBER, PLOTMETHOD, INVERTX,INVERTY, XMIN, XMAX, 2YMIN, YMAX, XTICSPACE, YTICSPACE, XBOTTOMLEFT, YBOTTOMLEFT ,XTOPRIGHT, 3YTOPRIGHT, PLOTZCOLOR, BORDERCOLOR, WINDOWCOLOR, AXISCOLOR, LABELCOLOR, 4PLOTCOLOR, TITLECOLOR, XAXISTITLE, YAXISTITLE, TITLE, TITLEONLY, 5FONTHEIGHT, FONTWIDTH, SCRONLY, Y2PLOT, Y3PLOT, PLOT3COLOR, X3, Y3, NUM3, 63W, PLOTZMETHOD) C WAIT FOR USER TO HIT RETURN C READ(*.*) CC CALL ENDGRAPHICS() OFLAG=0 RETURN ENDN OF MODULE SUBROUTINE OPENING C####################################################################### OOOOOOJFX'IPII’fiiI-fl-X'fl-X'I'I-I’ 000 w:~*n-**=+*x~*x-** 10 SUBROUTINE PLOTXY ******************************** ******************************** *x-sw:9s1~*s~*n~** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SUBROUTINE PLOTXY (TIME,FX,FY,Fz,PR,FXMAx,FYMAx,FZMAx,NP,Tx,TY,Tz, 1WEIGHT,TWEIGHT,CFALG,K,NF,FXMIN,FYMIN,FZMIN,TXM,TYM,N,FNAME, 2EWEIGHT,BW) DIMENSION TIME(512), FX(512), FY(512), FZ(512), PR(512) INTEGER WEIGHT(5) INCLUDE 'PLOT.DIF' LOGICAL*1 INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lPLOTMETHOD,INVERTX,INVERTY,TITLEONLY,SCRONLY,Y2PLOT,y3plot, 2PLOT2METHOD INTEGER*2 BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 1TITLECOLOR,NUMBER,PLOT2COLOR,FONTHEIGHT,FONTWIDTH,PLOT3COLOR CHARACTER TITLE*280,XAXISTITLE*80,YAXISTITLE*80,FNAME1*3,FNAME2*3, 1FNAME3*3,FNAME4*3,FNAME*12,TITLE1*80,TITLE2*80,TITLE3*80,TITLE4* 280,TITLE5*80,TITLE6*80,TITLE7*80,TITLE8*80,TITLE9*80,TITLE10*80 NUMBER=N INITIALISE THINGS FOR THE FIRST PLOT INITGRAPHICS=.TRUE. INITBORDER=.TRUE. INITWINDOW=.TRUE. INITAXES=.TRUE. INITFONTS=.TRUE. PLOTMETHOD=.TRUE. PLOT2METHOD=.TRUE. TITLEONLY=.FALSE. SCRONLY=.FALSE. Y2PLOT=.TRUE. FONTHEIGHT=10 FONTWIDTH=5 SET UP THE FANCY COLORS DEPENDING ON THE VALUE OF CFLAG (0 OR 1) BORDERCOLOR=14 IF (CFLAG.EQ.0.0) THEN WINDOWCOLOR=8 CFLAG=1.0 ELSE WINDOWCOLOR=1 CFLAG=0.0 END IF AXISCOLOR=10 LABELCOLOR=11 PLOTCOLOR=14 PLOT2COLOR=4 PLOT3COLOR=4 TITLECOLOR=15 IF(BW.EQ.1.0) THEN BORDERCOLOR=3 WINDOWCOLOR=15 AXISCOLOR=0 LABELCOLOR=0 PLOTCOLOR=0 PLOT2COLOR=1 PLOT3COLOR=3 TITLECOLOR=4 END IF Set file name Component READ (FNAME,10) FNAME1,FNAME2,FNAME3,FNAME4 FORMAT (A2,A1,A1,1x,A3) 000 000 0000000000 000 20 0000 0000000000 000 000 122 ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. FORCE IN THE x-DIR ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SET UP THE PLOT LIMITS for the screen part only PRED=2.0 XMIN=1.0 XMAX=6.0 XTICSPACE=2.0 Call scale subroutine to set Ymin, Ymax, Yticspace CALL SCALE (FXMIN,FXMAX,YMIN,YMAX,YTICSPACE) POSITION THIS PLOT AT THE TOP LEFT OF THE SCREEN XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.50 XTOPRIGHT=0.50 YTOPRIGHT=1.00 GIVE IT.A TITLE TITLE='FX AS A RATIO TO WEIGHT ' XAXISTITLE=' T I M E ( s e C . )' YAXISTITLE='FORCE RATIO' PLOT THE SECOND SET OF DATA (THE PROMPT) IF ((FNAME2.EQ.'t').OR.(FNAME2.EQ.'T')) THEN PR(NP)=FXMAX/PRED ELSE DO 20 I=1,N IF (PR(I).GT.0.0) PR(I)=FXMAX/PRED CONTINUE END IF READY TO CALL PROGRAM PLOT.FOR CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lTIME,FX,PR,NUMBER,PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAX, 2XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,X3,Y3,NUM3,BW 6,PLOT2METHOD) ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 2. FORCE IN THE Y-DIR ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SET UP THE PLOT LIMITS XMIN=1.0 XMAX=6.0 XTICSPACE=2.0 Call scale subroutine to set Ymin, Ymax, Yticspace CALL SCALE (FYMIN,FYMAX,YMIN,YMAX,YTICSPACE) POSITION THIS PLOT AT THE TOP RIGHT OF THE SCREEN XBOTTOMLEFT=0.50 YBOTTOMLEFT=0.50 XTOPRIGHT=1.00 YTOPRIGHT=1.00 GIVE IT A.TITLE ' TITLE='FY AS.A RATIO TO WEIGHT ' XAXISTITLE=' T I M E ( s e c . ) 0000 000000 0000 0000 0 000 000 0000 0000 30 40 123 YAXISTITLE='FORCE RATIO' Prevent a new screen INITGRAPHICS=.FALSE. INITFONTS=.FALSE. INITAXES=.TRUE. PLOT THE SECOND SET OF DATA (THE PROMPT) IF ((FNAME2.EQ.'t').OR.(FNAME2.EQ.'T')) THEN PR(NP)=FYMAx/PRED ELSE DO 30 I=1,N IF (PR(I).GT.0.0) PR(I)=FYMAx/PRED CONTINUE END IF READY TO CALL PROGRAM PLOT.FOR CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, 1TIME,FY,PR,NUMBER,PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAX, 2XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,X3,Y3,NUM3,BW 6,PLOT2METHOD) ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 3. FORCE IN THE z-DIR ++++++++++++++++ ++++++++++++++++++++++++++++++++++ Call scale subroutine to set Ymin, Ymax, Yticspace CALL SCALE (FZMIN,FZMAX,YMIN,YMAX,YTICSPACE) POSITION THIS PLOT AT THE BOTTOM LEFT OF THE SCREEN XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.00 XTOPRIGHT=0.50 YTOPRIGHT=0.50 GIVE IT A.TITLE TITLE='FZ AS.A RATIO TO WEIGHT XAXISTITLE=' T I M E ( s e c . )' YAXISTITLE='FORCE RATIO' PLOT THE SECOND SET OF DATA (THE PROMPT) IF ((FNAME2.EQ.'t').OR.(FNAME2.EQ.'T')) THEN PR(NP)=FZMAX/PRED ELSE DO 40 I=1,N IF (PR(I).GT.0.0) PR(I)=FZMAX/PRED CONTINUE END IF READY TO CALL PROGRAM PLOT.FOR INITAXES=.TRUE. CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, R NUMBER PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAx, 2§T¥§SPAC§,YTICSPACE,XBOTTOMLEFT,YEOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,X3,Y3,NUM3,EW 6,PLOT2METHOD) 0000000 000 0 00000 OOOOOOfi-X'Ififlki-fl'li'i-‘f’ffl'ifoo 50 60 70 80 90 100 110 130 140 # ** * * **=+*l-**Ik*I-*i-W# ### 124 ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 4. OUTPUT VALUES ++++++++++++++++ ++++++++++++++++++++++++++++++++++ just plot the title on this one XBOTTOMLEFT=0.50 YBOTTOMLEFT=0.00 XTOPRIGHT=1.00 YTOPRIGHT=0.50 TITLE=' FORCE—RATIOS ' WRITE (TITLE1,50) FNAME,FNAME3 FORMAT ('FILE = ',A12,2x,'WITH ',A2,2X,'PER.') IF ((FNAME2.EQ.'t').OR.(FNAMEZ.EQ.'T')) THEN WRITE (TITLE3,60) FNAME1,FNAME4 ELSE WRITE (TITLE3,70) FNAME1,FNAME4 END IF FORMAT ('TRANSIENT MOTION = ',A3,2x,'TEST = ',A3/) FORMAT ('PERIODIC MOTION = ',A3,2x,'TEST = ',A3/) WRITE (TITLE4,80) (WEIGHT(I),I=1,5) FORMAT (”Weights = ',5(I3,',')) WRITE (TITLE5,90) TWEIGHT,EWEIGHT FORMAT ('Tota1= ',F4.0,' ESTIMATE= ',F7.2,' lbs') Find the max. value of force in the x,y directions (absolute force IF (ABS(FXMIN).GT.ABS(FXMAX)) FXMAX=FXMIN IF (ABS (FYMIN) .GT.ABS (FYMAX) ) FYMAX=FYMIN IF (ABS(FXMIN).GT-ABS(FXMAX)) TX=TXM IF (ABS(FYMIN).GT.ABS(FYMAX)) TY=TYM WRITE (TITLE6,100) K,NF WRITE (TITLE7,110) FXMAX,TX WRITE (TITLE8,120) FYMAX,TY WRITE (TITLE9,130) FZMAx,Tz WRITE (TITLE10,140) FORMAT (IZ,' of '13,2x,' FMAx TIME') FORMAT ('x—dir',3x,F8.3,6x,F6.3) FORMAT ('Y-dir',3X,F8.3,6X,F6.3) FORMAT ('z-dir',3x,F8.3,6x,F6.3) FORMAT ('HIT ENTER TO CONTINUE....’) CALL OUTPUTS (XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT,TITLE, lTITLEl,TITLE2,TITLE3,TITLE4,TITLE5,TITLE6,TITLE7,TITLE8,TITLE9, 2TITLE10,CFLAG,BW) WAIT FOR USER TO HIT RETURN ( only for testing the program , otherwise let the program run With no stops READ (*,*) IF (K.EQ.NF) CALL ENDGRAPHICS () RETURN ENEND OF SUBROUTINE PLOTXY ### at#WE1*TEE#FEEWWEEEWEEEETEEEEEEE . 3(- a» 2(- s * :(- I? :(~ 1» 1' * If If at u- if at u» w l- * n- x- * l- * up 2(- a; s :1- 2(- :- SUBROUTINE MAX ***************** *************** *************** **************~k** MAY 05,1992 FINDING THE MAX. VALUE OF DATA ** **-**-*t-**-* ++++++++++++++++ 1. SUBROUTINE MAX (DATA,FMAX,T,TIMESTEP) 000 000 OOOQOO*************OO OEEEEEEEEEEEEOOO 10 1 am: *l-fl'fl-II’X'iffl-fl-X'I-I'fitfit f * * * * * * * * * * * * O #43: # #£################################################################### 125 TO FIND THE MAX. FORCE DIMENSION DATA(512) T=0.0 N=0 =-1E14 NMAX=0 DO 10 I=1,512 N=N+1 IF (DATA(I).GT.FMAX) THEN FMAX=DATA(I) NMAX=N GO TO 10 END IF CONTINUE T=REAL(NMAX)*TIMESTEP RETURN ENDN OF SUBROUTINE * i. E E * E E E E E E f # E E E * E E i I" E E E E E E E E E E E E E E E * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ***************************************************************** SUBROUTINE MIN MAY 05,1992 FINDING THE MIN. VALUE OF DATA Ex-Ex-EEarEx-EE ++++++++++++++++++++++++++++++++++ 1. DATA.INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ ++++++++++++++++ SUBROUTINE MIN (DATA,FMIN,T,TIMESTEP) TO FIND THE MIN. FORCE DIMENSION DATA(512) T=0.0 N=O FMIN=0.0 NMAX=0 DO 10 I=1,512 N=N+1 IF (DATA(I).LT.FMIN) THEN FMIN=DATA(I) NMAX=N GO TO 10 END IF CONTINUE T=REAL(NMAX)*TIMESTEP RETURN END UBROUTINE MIN END OF S #################################################################### f#f*f# .#f#f*f”f”f#f*f#f#f#f#f*f#f#2#f#f#f#f#f#f#f*f#f*f#f#f*f#f*f#f# * * * * * * * * * * * * * * * * * * * * * * * * * * * k * * * * * * * SUBROUTINE outputs : * * * k * * * * * * * * * * * * * * * * * i * * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * MAY 15,1992 Plotting the final values of Forces : * * * * * * * * * * * k * * * * * * * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * START OF MODULE outputs.for INCLUDE 'FGRAPH.FI' 0 00000 0000 000 0000 10 20 126 SUBROUTINE OUTPUTS (XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 1TITLE,TITLE1,TITLE2,TITLE3,TITLE4,TITLE5,TITLE6,TITLE7,TITLE8, 2TITLE9,TITLE10,CFLAG,BW) INCLUDE 'FGRAPH.FD' INTEGER*2 DUMMY,XWIDTH,YHEIGHT,IX1,IY1,IX2,IY2,IX,IY,FONTHEIGHT, 1FONTWIDTH REAL*4 XWl,YW1,XW2,YW2,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT, 1YTOPRIGHT,YT CHARACTER TITLE*80,FONTCOMMAND*10,FONTSTRING*7,FONTFILE*12,TITLE1* 180,TITLE2*80,TITLE3*80,TITLE4*80,TITLE5*80,TITLE6*80,TITLE7*80, 2TITLE8*80,TITLE9*80,TITLE10*80 RECORD / VIDEOCONFIG / SCREEN RECORD / WXYCOORD / WXY RECORD / XYCOORD / POSITION COMMON SCREEN SET UP FONT FOR USE FONTFILE='HELVB.FON' FONTCOMMAND="T'HELV'" IF (REGISTERFONTS(FONTFILE).LT.0) THEN WRITE (*,10) FONTFILE FORMAT (' FONT FILE ',A12,' NOT FOUND IN PLOT') STOP END IF FONTHEIGHT=28 FONTWIDTH=16 WRITE (FONTSTRING,20) FONTHEIGHT,FONTWIDTH FORMAT ('H',I2.2,'W",Iz.2,'B') DUMMY=SETFONT(FONTCOMMAND//FONTSTRING) GET SCREEN DIMENSIONS XWIDTH=SCREEN.NUMXPIXELS YHEIGHT=SCREEN.NUMYPIXELS SET A VIEW PORT, WITH CORNER COORDINATES DEFINED As A FRACTION OF THE TOTAL SCREEN SIZE (0.0,o.0) IS BOTTOM LEFT OF SCREEN (1.0,1.0) IS TOP RIGHT OF SCREEN TRANSLATE CORNER COORDINATES TO PIXEL COORDINATES IX1=INT(XBOTTOMLEFT*FLOAT(XWIDTH)) IX2=INT(XTOPRIGHT*FLOAT(XWIDTH)) IY1=YHEIGHT—INT(YTOPRIGHT*FLOAT(YHEIGHT)) IY2=YHEIGHT—INT(YBOTTOMLEFT*FLOAT(YHEIGHT)) CALL SETVIEWPORT (IX1,IY1,IX2,IY2) CREATE.A REAL COORDINATE WINDOW XW1=0.0 YW1=0 XW2=1.0 YW2=1.0 DUMMY=SETWINDOW(.TRUE.,XW1,YW1,XW2,YW2) SET THE COLOR OF THE BACKGROUND WINDOW DEPENDING ON THE FLAG "CFLAG" ... NOTE THE SETUP HAS TO BE OPPOSITE TO THE PLOTXY IF (CFLAG.EQ.1.0) THEN DUMMY=SETCOLOR(8) ELSE DUMMY=SETCOLOR(1) END IF IF(BW.EQ.1.0)DUMMY=SETCOLOR(15) DUMMY=RECTANGLE W($GFILLINTERIOR,XW1,YW1,XW2,YW2) DRAW.A BORDERIAROUND THE VIEWPORT DUMMY=SETIOU?RUMMY SETCOLOR(3) IF BW.E . . D = DUMMY=RECTANGLE;W($GBORDER,XWl,YW1,XW2,YW2) This is the screen design of the data output TITLELENGTH=50 POSITION THE TITLE CENTERED (MORE OR LESS) ABOVE GRAPH YT=0.9*(YW2+YW1) CALL MOVETO W (0.9*(XW2+XW1),YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX—(TITLELENGTH*FONTWIDTH)/2 IY=IY-FONTHEIGHT/2 000 00000 000M 0000 127 enter the program name (Title: from the main program) DUMMY=SETCOLOR(15) CALL MOVETO (IX,IY,POSITION) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY+2,POSITION) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) CHANGE THE FONT TO.A FIXED FONT TO CREATE THE FORCES' TABLE HERE, COURIER FONTS WILL BE USED. FONTFILE='courb.fon' FONTCOMMAND="T'courier'" FONTHEIGHT=12 FONTWIDTH=9 IF (REGISTERFONTS(FONTFILE).LT.0) THEN WRITE (*,10) FONTFILE STOP END IF WRITE (FONTSTRING,20) FONTHEIGHT,FONTWIDTH DUMMY=SETFONT(FONTCOMMAND//FONTSTRING) YT=0.8*(YW2+YW1) CALL MOVETO W (0.73*(XW2+XW1),YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-(TITLELENGTH*FONTWIDTH)/2 IY=IY—FONTHEIGHT/2 DUMMY=SETCOLOR(0) CALL MOVETO (IX,IY+10,POSITION) CALL OUTGTEXT (TITLE1(1:TITLELENGTH)) CALL MOVETO (IX,IY+25,POSITION) CALL OUTGTEXT (TITLE2(1:TITLELENGTH)) CALL MOVETO (IX,IY+40,POSITION) CALL OUTGTEXT (TITLE3(1:TITLELENGTH)) CALL MOVETO (IX,IY+60,POSITION) CALL OUTGTEXT (TITLE4(1:TITLELENGTH)) CALL MOVETO (IX,IY+75,POSITION) CALL OUTGTEXT (TITLE5(1:TITLELENGTH)) IF(BW.EQ.1.0)GOTO 22 DUMMY=SETCOLOR(11) CALL MOVETO (IX,IY+11,POSITION) CALL OUTGTEXT (TITLE1(1:TITLELENGTH)) CALL MOVETO (IX,IY+26,POSITION) CALL OUTGTEXT (TITLE2(1:TITLELENGTH)) CALL MOVETO (IX,IY+41,POSITION) CALL OUTGTEXT (TITLE3(1:TITLELENGTH)) DUMMY=SETCOLOR(10) CALL MOVETO (IX,IY+61,POSITION) CALL OUTGTEXT (TITLE4(1:TITLELENGTH)) CALL MOVETO (IX,IY+76,POSITION) CALL OUTGTEXT (TITLE5(1:TITLELENGTH)) CONTINUE Change font size FONTHEIGHT=10 FONTWIDTH=8 WRITE (FONTSTRING,20) FONTHEIGHT,FONTWIDTH DUMMY=SETFONT(FONTCOMMAND//FONTSTRING) Filling the Table of Forces DUMMY=SETCOLOR(14) IF(BW.EQ.1.0)DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY+110,POSITION) CALL OUTGTEXT (TITLE6(1:TITLELENGTH)) DUMMY=SETCOLOR(15) IF(BW.EQ.1.0)DUMMY=SETCOLOR(0) CALL MOVETO (IX,IY+130,POSITION) CALL OUTGTEXT (TITLE7(1:TITLELENGTH)) 000 **I+*I-**-**-**-** O()O()*i-**l$*i-*i'**-** 0000000 * * 128 CALL MOVETO (IX,IY+145,POSITION) CALL OUTGTEXT (TITLE8(1:TITLELENGTH)) CALL MOVETO (IX,IY+160,POSITION) CALL OUTGTEXT (TITLE9(1:TITLELENGTH)) DUMMY=SETCOLOR(0) CALL MOVETO (IX,IY+180,POSITION) CALL OUTGTEXT (TITLE10(1:TITLELENGTH)) DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY+181,POSITION) CALL OUTGTEXT (TITLE10(1:TITLELENGTH)) ALL DONE CALL UNREGISTERFONTS () RETURN END *********** ** ******** ********* ******** SUBROUTINE SCALE ******************* ******************* AUG 22,1992 * * ** ** VERTICAL PLOT SCALE ‘1:******************************** ********************************* ***************************************************************** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ ++++++++++++++++++++++++++++++++++ ++++++++++++++++ SUBROUTINE SCALE (FMIN,FMAX,ST,FIN,TICK) REAL*8 U,C TT=ALOG10(FMAx-FMIN) MT=IRINT(TT) TT=10.**(TT-MT) IF (TT.LE.2) THEN TT=2. ELSE IF (TT.GT.5) THEN TT=10. ELSE TT=5. END IF TICK=TT*10.**(MT—l) U=1.00001*TICK C=DLOG10(U) IF (C.LT.0) THEN L=IDINT(C)~1 ELSE L=IDINT(C) END IF C=DBLE(10.)**L TICK=C*IDNINT U/C) N1=IRINT(FMIN TICK) N2=IRINT(FMAX/TICK)+1 ST=N1*TICK FIN=N2*TICK RETURN 1. DATA INPUT/OUTPUT FIN = Ending value; TICK 1-##-** **-*1-** EE -----—-—-—-—-----------—------—--------------—------—_----—-—-—-—--—-~_ This routine determines the starting and ending values to be used in a full scale plot. OUTPUT : ST = Starting value; THIS IS A.MODIFIED VERSION OF.A PROGRAM WRITTEN ORIGINALLY BY R. HARIC ——_-————-—-—-————-.————-u-————————————u——————---—~————-————_——_—————-c———-- 129 $ SRECTRA.FOR STORAGE: SUBROUTINE SPECTRA * * * * * * * * * * * * * * * * * * * * * * * * * * * k * * * * * AUG 22,1992 AUTOCORRELATION AND AUTOSPECTRA PLOTS * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * k * * i * * * k k * * * * * * * k * * * * * * * * * * * * *I-*I-**-**>+*i-**#H# f * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * k 'k * * * * ***************************************************************** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ ——_————-—..---‘_-———————-----—-——_———————-———--——--_--——————————————-—_—_ Description: PROGRAM TO ESTIMATE THE AUTOCORRELATION, AUTOSPECTRA, AVERAGE AUTOCORRELATIONC, AND AVERAGE AUTOSECTRUM OF FORCE-TIME SERIES. IT IS A.MODIFICATION OF A PROGRAM ORIGINALLY WRITTEN BY: Ronald S. Harichandran, and E. Vanmarke, "Space-Time Variation of Earthquake Ground Motion", Research Report R84—12, Dept. of Civil Engineering, MIT, 1984. ————~--———————--——--—————-—-—-——-——-——-—-——-————-—————-—--—————————-——— The autocorrelation function and the one—sided, normalized (i.e. area) autospectrum of each force-series are stored in files with the same name as the force-series and with filetype .ACN and .AS respectively. The average of all the autocorrelation functions a of the normalized autospectra are stored in files named AVERAGE. and AVERAGE.ASP; If autospectra are estimated, then a file name SUMMARY.OUT, containing useful summary tables, is created. Description of input: Input are received from a file named "PARAM.DAT" which contains: (1) General Informatio about the run: N No. of Samples FS Sampling Rate (2) For autospectral estimation: NA Number of force-series Name(L,I) Names of force-series (Each in a separate lin (3) Input for autospectral computations: M Section size (must be a power of 2), 2 <= M <= 2 IWIN Window t e, l for rect. window, 2 for Hamming wi L Half-widt of lag window, 2 <= L<=(M/2+1); NFFT FFT size used for the spectral estimate, (2L-l) OF Output frequency (<= fs/2) C * * * 1|: * 1k * * ~k * * 1k * C C C C C C C C c C C C C C c C E Description of output: c C c C C C C C C c C C C C C C C C c C C C C _______________________________________________________________________ SUBROUTINE SPECTRA (LL,K,NA,NAME,OF,RMS,XMEAN,BW,noavg) COMPLEX SXY(513) REAL SACORR(1025),SASPEC(513) REAL SX(513),CCF(1025),XA(1025),AXIS(1025),NULL(1) CHARACTER FILE1*12,C*80,XLAB*80,YLAB*80,INPUTS*40,0UTPUTS*40 CHARACTER*8 NAME(20,200) COMMON /FILEs/ INPUTS,OUTPUTS,LINP,LOUT COMMON /CMP/ M,N,FS,IWIN,L,NFFT,DF COMMON /CMPD1/ MAXM COMMON /OUT1/ C COMMON /SM/ FMAX,NP DATA (SACORR(I),I=1,512)/512*0./,(SACORR(I),I=513,1025)/513*0./ DATA.(AXIS(I),I=1,1025)/1025*0./ DATA.SASPEC/513*0./ C C Begin spectral estimation 000 000 10 000 20 30 CCC CCC 4O CCC ccc 000 50 C 130 DF=Fs/FLOAT(NFFT) NF=NFFT/2+1 XA(1)=0. XA(2)=FMAX- DF CALL AUTOSCL (XA, 2, FTICK, FST, FMAX) NFMAX= FMAX/DF+1. 00001 NOF=ANINT(OF/DF)+1 IF (NOF.GT.NF) NOF=NF Compute the autocorrelation functions and autospectra. II=K OPEN (6,STATUS='OLD',FILE=OUTPUTS(1:LOUT)//NAME(LL,II)//'.FRC') skip information lines at the top of the input file (5 lines) READ (6,10) FORMAT (5(/)) CALL CMPSE (CCF,SXY,VAR,PA,XMEAN) CLOSE (6) VAR=SQRT(VAR) RMS=VAR ADD THE MEAN TO THE PEAK (IT WAS REMOVED IN THE SPECTRUM CALC.) IF (PA.GT.0) THEN PA:PA+XMEAN ELSE PA=PA—XMEAN END IF DO 20 J=1, M/2+1 XA(J)=J-1 SACORR(J)=SACORR(J)+CCF(J)/NA WRITE (FILE1,30) NAME(LL,II),'.ACN' FORMAT (8A) ='AUTOCORRELATION ' OPEN (8, FILE=OUTPUTS(1: LOUT)//FILE1) CALL OUT (8, M/2+1, CCF, 0.,1.) XLAB= 'Lag' YLAB='Autocorrelation' XTICK=0. YTICK=0. CALL AUTOSCL (XA, M/2+1, XTICK, XST, XFIN) CALL AUTOSCL (CCF, M/2+i, YTICK, YST, YFIN) CALL PLTS (0.,FLOAT(M/2+1), YST, 1.,XTICK, YTICK, XLAB, YLAB, C, AXIS, XA, 1CCF, M/2+1,11,BW) WRITE (FILE1, 30) NAME(LL, II),' .ASP' ='AUTOSPECTRUM DO 40 J=1, NF SX(J)=CABS(SXY(J)) XA(J)=(J- 1)*D F SASPEC(J)=SASPEC(J)+SX(J)/NA OPEN (9, FILE= -OUTPUTS(1: LOUT)//FILE1) CALL OUT (9, NOF, sx, 0.,DF) XLAB='Frequency (Hz)' YLAB='Spectral Density' YTICK = -2. YTICK=0. CALL AUTOSCL (SX, NFMAX, YTICK, YST, YFIN) CALL PLTS (0.,FMAX, YST, YFIN, FTICK, YTICK, XLAB, YLAB, C, AXIS, XA, SX, 1NFMAX, 2, BW) call plts for the final information on the screen =' S P E c T R A.L A.N A L Y s I 5' CALL PLTS (0.,1. o.,1.,1 1. ,XLAB, YLAB, C, AXIS, NULL, NULL, 1, 5, BW) WRITE (C, 50) NAME(LL, II), RMS, PA, XMEAN lFORMAT (1x, FILE: ,A8, T25, 'RMS:',F9.4,T42,'PEAK:',F7.2,T62,'MEA N: F7.2) CALL PLTS (0.,1.,0.,1.,l.,1.,XLAB,YLAB,C,AXIS,NULL,NULL,1,7,BW) g Output average correlations and spectrum if(noavg.eq.0) then 60 70 000 131 IF (II.EQ.NA) THEN FILE1='ANERAGE. ACN' C='AVERAGE CORRELATION OPEN (10, FILE=OUTPUTS(1: LOUT)//FILEl) CALL OUT (10, M/2+1, SACORR, 0.,1.) DO 60 J=1, M/2+1 XA(J)=J-1 XLAB='Lag' YLAB='Autocorrelation' XTICK=O. YTICK=O. CALL AUTOSCL (XA.,M/2+1, XTICK, XST, XFIN) CALL AUTOSCL (SACORR, M/2+1, YTICK, YST, YFIN) CALL PLTS (O. ,FLOAT(M/2+l),YST,1.,XTICK,YTICK,XLAB,YLAB,C,AXIS,XA, lSACORR, M/2+1,3 3,BW) FILE1='AVERAGE.ASP' C='AVERAGE SPECTRmM ' OPEN (11, FILE=OUTPUTS(1: LOUT)//FILE1) CALL OUT (11, NF, SASPEC, 0. ,DF) DO 70 J=1, M/2+1 XA(J)= (J—1)*DF XLAB='Frequency (Hz)' YLAB='Spectra1 Density' YTICK=O. CALL AUTOSCL (SASPEC,NFMAX,YTICK,YST,YFIN) CALL PLTS (0.,FMAX,YST,YFIN,FTICK,YTICK,XLAB,YLAB,C,AXIS,XA, 13ASPEC,NFMAX,4,BW) call plts for the final information on the screen =' A v E R.A G E A.N A L Y S I 3' CALL PLTS (0.,1.,0 .,1.,1.,1. ,XLAB, YLAB, c, AXIS, NULL, NULL, 1, 6, BW) =' E N D O F F I L E P R O C E S s I N G' CALL PLTS (0.,1.,0.,1.,1.,1.,XLAB, YLAB, c,AXIs, NULL,NULL,1,8,BW) END IF END IF END SUBROUTINE OUT (IF,NN,X,XMIN,XINC) ——————--_-———n—-_——-———————-——-———-———————_——_————————_—-—————————_— C This subroutine outputs an array to a file. I C____. 10 20 30 4O 00 #41: 00000000 ##3##: t1! —-_.. ~-———_-————-—-——---————-—--—————————-————_————-—--—-—-—————————————— REAL X(*) CHARACTER*80 c COMMON /CMP/ M,N,FS,IWIN,L,NFFT,DF COMMON /OUT1/ c FORMAT (A80/) FORMAT ' M= ,I5, '; No. of Samples =' ,IS, ; Sampling Frequency éf‘,g?.l ' Window = Hamming' ,'; ‘Window half-width =' I4, ; NFFT = ,I FORMAT ' M =' ,15, '; NO. of Sm mples =' ,IS, Sampling Frequency 1=',F6.l ' Window = Rectangular' , Window half-width = ,I4, NF 2FT =',I5) OPEN (IF,STATUS='UNKNOWN') WRITE (IF,10) 0 IF (IWIN.EQ.1) WRITE (IF,30) M,N,FS,L,NFFT IF (IWIN.EQ.2) WRITE (IF,20) M,N,FS,L,NFFT WRITE (IF,*) DO 40 I=1,NN WRITE (IF *) XMIN+(I-1)*XINC,X(I) CLOSE (IFS RETURN ND ####################################### ####################################### SUBROUTINE CMPSE (CCF,SXY,XVAR,PA,XMEAN ~w¢#= -———————_--—-——-_-———-—--————————-—-_———————--——-—————-——————q————-—-- This subroutine estimates the auto correlation functions, and the normalized auto amplitude spectra. The correlation method for power spectrum estimation is used. Thi modification of the program written by: L. Rabiner, Bell Laboratories, MUrray Hill, New Jersey 07974 R. Schafer, et. a1, Georgia Institute of Technology, Atlanta, Georgi 132 This method is based on the technique described by C. M. Rader, In the IEEE Trans. on Audio and Elect., VO1. 18, No. 4, pp 439-442 Input: M is the section size (must be a power of 2); 2 = M <= 1024. N is the number of samples to be used in the analysis. MODE is the data format type; _MDDE = 0 for autocorrelation and autospectrum; FS 18 the sampling frequency in Hz (i.e., number of samp second). IWIN is the window type; IWIN = 1 for rectangular window; IWIN = 2 for Hamming window. (2L-1) is the window base width used in the spectral est 2 <= L <= (M/2+1). NFFT is the FFT size used in the spectral estimate; (2L-1) <= NFFT <= 1024. Variables returned to calling program: CCF(M/2+1) = Autocorrelation function for +ve lags; SXY(NFFT/2+1) = Autospectrum (take absolute of SXY); XVAR = variance of the force-series; PA = Absolute peak of the force-series. REAL XA(1025),CCF(*) COMPLEX X(1025),Z(1025),XMN,XI,YI,SXY(*) COMMON /CMP/ M,N,FS,IWIN,L,NFFT,DF PARAMETER (PI=3.141592654) 000OOOOOOOOOOOOOOOOOOOOOO c C DEFINE CONSTANTS. NSECT IS THE TOTAL NUMBER OF ANALYSIS SECTIONS. c LSHFT IS THE SHIFT BETWEEN ADJACENT ANALYSIS SAMPLES C LSHFT=M/2 MHLF1=LSHFT+1 NSECT=(FLOAT(N)+FLOAT(LSHFT)-l.)/FLOAT(LSHFT) c C 33 IS THE GENERATOR SAMPLE NUMBER. AT THE VERY FIRST CALL TO SUBROUTI c GETX THIS IS SET AT zERo SO THAT ANY NECESSARY INITIALIZATION C CAN BE PERFORMED IN THESE SUBROUTINES. c NRD IS NUMBER OF SAMPLES OF GENERATOR OUTPUT To BE COMPUTED OR READ. c SS=0. NRD=LSHFT XSUM=0. YSUM:0. XXSUM?0. YYSUM=0. PA=0. LOOP TO CALCULATE MEANS AND VARIANCES OF X AND Y DATA USE GETX TO READ NRD SAMPLES FROM X GENERATOR STARTING AT SAMPLE SS 0000 DO 30 K=1,NSECT IF (K.EQ.NSECT) NRD=N-(K-1)*NRD CALL GETX (XA,NRD,SS,N) DO 10 I=1,NRD XSUM:XSUM+XA(I) 10 XXSUM?XXSUM+XA(I)*XA(I) C C COMPUTE THE ABSOLUTE MAXIMUM OF THE NRD VALUES OF XA. C DO 20 I=1,NRD 20 IE‘ (ABS(XA(I)).GT.ABS(PA)) PA=XA(I) IF (K.EQ.1) SS=1. SS=SS+ELOAT(NRD) 30 CONTINUE FN=FLQAT(N) XMEAN=XSUM/FN XVAR=XXSUM/FN-XMEAN*XMEAN SQVAR=SQRT(XVAR*XVAR) XMN=CMPLX(XMEAN,XMEAN) REMOVE THE MEAN FROM THE PEAK (MEAN WILL BE A DELTA FUNCTION AT ZE 000 IF (PA.GT.0) THEN PAFPA-XMEAN 133 ELSE PA=PA+XMEAN END IF C C LOOP TO ACCUMULATE CORRELATIONS C SS=1. NRDY=M NRDX=LSHFT DO 40 I=1,NHLE1 40 Z(I)=(0.,O.) DO 110 K=1,NSECT NSECT1=NSECT-l IF (K.LT.NSECT1) GO TO 60 NRDY=N—(K-l)*LSHFT IF (K.EQ.NSECT) NRDx=NRDY IF (NRDY.EQ.M) Go TO 60 NRDY1=NRDY+1 DO 50 I=NRDY1,M 50 X(I)=(0.,0.) c c READ NRDY SAMPLES FROM x GENERATOR STARTING AT SAMPLE SS C 60 CALL GETx (XA,NRDY,SS,N) DO 70 I=1,NRDY 70 X(I)=CMPLX(XA(I),XA(I)) DO 80 I=1,NRDY 80 X(I)=X(I)-XMN NRDx1=NRDx+1 DO 90 I=NRDX1,M 9o X(I)=CMPLX(0.,AIMAG(X(I))) c c CORRELATE x AND Y SECTIONS C DO EVEN-ODD SEPARATION AND ACCUMULATE CONJG(X)*Y C CALL FFT (X,M,O) DO 100 I=2,LSHFT J=M+2—I XI=(X(I)+CONJG(X(J)))*.5 YI=(X(J)-CONJG(X(I)))*.5 YI=CMPLX(AIMAG(YI),REAL(YI)) Z(I)=Z(I)+CONJG(XI)*YI 100 CONTINUE XI=X(1) Z(1)=Z(1)+CMPLX(REAL(XI)*AIMAG(XI),0.) XI=X(MHLF1) Z(MHLF1)=Z(MHLF1)+CMPLX(REAL(XI)*AIMAG(XI),0.) SS=SS+FLOAT(LSHFT) 110 CONTINUE C c INVERSE DFT TO GIVE CORRELATION C DO 120 I=2,LSHFT J=M+2—I X(I)=Z(I) X(J)=CONJG(Z(I)) 120 CONTINUE X(1)=Z (l) X(MHLF1)=Z(MHLF1) C CALL FFT (X,M,1) g STORE NORMALIZED CORRELATIONS FOR POSITIVE LAGS. DO 130 I=1,MHLE1 XA(I)=REAL(X(I))/FN CCF(I)=XA(I)*M/SQVAR XA(I)=CCF(I) C130 CONTINUE C WINDOW CORRELATION USING (2L-l) POINT WINDOW. FOR AUTOCORRELATION USE C OF SYMMETRY. g In section below, K = —1 for negative lags and K = l for positive lags K=-1 J2=1 L2=L NLAST=NFFT~L2+1 134 IF (IWIN.EQ.2) THEN DO 140 I=J2,L2 140 éfigl§FXA(I)*(0.54+0.46*COS(PI*FLOAT(K*(I-l))/FLOAT(L-1))) C Put negative lag values at right end of XA(I) C IF (K.EQ.-l) THEN DO 150 I=2,L2 J=NFFT+2-I 150 XA(J)=XA(I) END IF C C DO 160 I=L2+1,NLAST 160 XA(I)=O. DO 170 I=1,NFFT 17o X(I)=CMPLX(XA(I),0.) C C COMPUTE NORMALIZED ONE-SIDED AUTOSPECTRUM OR CROSS SPECTRUM C CALL FFT (X,NFFT,0) DO 180 I=1,NFFT/2+l 180 SXY(I)=2.*X(I)/DF RETURN END C####################################################################### C####################################################################### SUBROUTINE AUTOSCL (X,N,TICK,ST,FIN) C This routine determdnes the starting and ending values to be used in P C for a full scale plot. ' C INPUT : X = X or Y array; N = Number of elements in X; C TICK = Distance between tick marks; a) >0 to use supplied valu C b) =0 to choose automatically C OUTPUT : ST = Starting value; FIN = Ending value; TICK (If originall —————-—u-——--——-_————-——-————————--—~-—‘--—-——————-—~————-----~-—-—-—_~ O FIN=X(1) ST=X(1) DO 10 I=1,N IF (X(I).GT.FIN) FIN=X(I) 10 IF (X(I).LT.ST) ST=X(I) IF (TICK.GT.0.) THEN ST=TICK*IRINT(ST/TICK) FIN=TICK*(IRINT(FIN/TICK)+1) ELSE IF (TICK.EQ.0) THEN TT=ALOGlO(FIN-ST) MT=IRINT(TT) TT=10.**(TT-MT) IF (TT.LE.2) THEN TT=2. ELSE IF (TT.GT.5) THEN TT=10. ELSE TT=5. END IF TICK=TT*10.**(MT-l) U=1.00001*TICK C=DLOG10(U) IF (C.LT.O) THEN L=IDINT(C)-l ELSE L=IDINT(C) END IF C=DBLE(10.)**L TICK=C*IDNINT(U/C) N1=IRINT(ST/TICK) N2=IRINT(FIN/TICK)+1 ST=N1*TICK FIN=N2*TICK NT=N2-Nl END IF 135 RETURN END INTEGER FUNCTION IRINT(X) IF (X.LT.0) THEN IRINT=INT(X)-1 ELSE IRINT=INT(X) END IF RETURN ################## ################## SUBROUTINE FFT (X, mm: #m¢ ##tm z:w# ~ #mt H3fl¢ Earn: ma: ~vma: am: am: *at C This subroutine computes the discrete Fourier Transform or the inverse C transform of X. C ——————————————————————————————————————————————————————————————————————— C Input: X = 2**M complex array that initially contains the input C and on return contains the transform; C N = 2**M points. (Must be power of 2, else infinite loop resu C INV = 0 for direct transform; 1 for inverse transform. ———————————--——-—-—-u—————_————-—--——_----—---——————————_———-———————_——— 00 COMPLEX X(*),U , MFALOG(FLOA T(N ) Nv2=N/2 NM1=N- 1 J=l DO 40 I=1,NM1 IF (I.GE.J) GO TO 10 T=X(J) X(J)=X(I) X(I)=T 10 K=NV2 20 IF (K.GE.J) GO To 30 J=J-K K=K/2 GO To 20 30 J=J+K 40 CONTINUE PI=4.0*ATAN(1.0) DO 70 L=1,M LE=2**L LE1= LE/z U=(1.0, 0.0) w=CMPLX(COS(PI/FLOAT(LEl)),-SIN(PI/ELOAT(LE1))) IF (INV. NE. 0) w=CONJG(W) DO 60 J=1,LE1 DO 50 I=J,N,LE IP=I+LE1 W,T )/ALOG(2.)+.1 T=X(IP) *U X(IP)=X(I)-T X(I)=X(I)+T 50 CONTINUE U=U*W 60 CONTINUE 7O CONTINUE IF (INV.EQ.1) RETURN DO 80 I= 1, N X(I)=X(I)/N 80 CONTINUE RETURN END ########################### ########################### SUBROUTINE GETX (X, NRD, SS, N #mt aw: ~«wsc _-_______—--_—____-—---_—__—_——----_--___——_—-_-_-_-.-—~__---_—--——_-_-— C This subroutine passes the required length of a TIME SERIES to the cal C program. On the first call to this subroutine SS = 0, and the ent C values of the force- series is read from logical unit 11, and the f C values are placed in array X. On all subsequent calls NRD values C force- series starting from sample SS are placed in array X. -——-—-~—---—-—_---——----‘—----_—————----_——--~——---——-—--—---——-——‘-~—- REAL X(*),XSTORE(7500) 10 FORMAT (45X,F13.6) OOOOOOO»**»»*******»OO 0000 000 N 0 :m# 136 J=SS IF (J.EQ.0) THEN §EID (6,10) (XSTORE(I),I=1,N) END IF DO 20 I=1,NRD X(I)=XSTORE(J+I-l) RETURN END . ###################################################################### ###################################################################### ******'k**************************** ********‘k************************** * * * SUBROUTINE PLTS * * * ***at-******************************* ******'k**************************** * 'k * AUG 22,1992 TO PLOT CORRELATION AND SPECTRUM : 'k * * * * * * * * * * * * i k * * * * * i * * * * * * * * * * * * * * * * k * * * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * ***************************************************************** * ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SUBROUTINE PLTS (XMIN,XMAX,YMIN,YMAX,XTICSPACE,YTICSPACE, lXAXISTITLE,YAXISTITLE,TITLE,AXIS,X,Y,NUMBER,NN,BW) LOGICAL*1 INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lPLOTMETHOD,INVERTX,INVERTY,TITLEONLY,SCRONLY,Y2PLOT,y3plot, 2PLOT2METHOD INTEGER*2 BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, lTITLECOLOR,NUMBER,PLOTZCOLOR,FONTHEIGHT,FONTWIDTH,PLOT3COLOR DIMENSION X(NUMBER), Y(NUMBER), AXIS(*) CHARACTER TITLE*280,XAXISTITLE*80,YAXISTITLE*80 INITIALISE THINGS FOR THE FIRST PLOT INITGRAPHICS=.TRUE. INITBORDER=.TRUE. INITWINDOw=.TRUE. INITAXES=.TRUE. INITFONTS=.TRUE. PLOTMETHOD=.TRUE. TITLEONLY=.FALSE. y2plot should be set to "true" to show the zero axis Y2PLOT=.TRUE. Y3PLOT=.FALSE. FONTHEIGHT=10 FONTWIDTH=5 SET UP THE FANCY COLORS BORDERCOLOR=14 WINDOWCOLOR=1 AXISCOLOR=10 LABELCOLOR=11 PLOTCOLOR=14 PLOT2COLOR=4 PLOT3COLOR=4 TITLECOLOR=15 IF(BW.EQ.1.0) THEN BORDERCOLOR=3 WINDOWCOLOR=15 AXISCOLOR=0 LABELCOLOR=4 PLOTCOLOR=0 PLOT2COLOR=1 PLOT3COLOR=3 TITLECOLOR=0 00000 00000 00000 00 00000 00000 137 END IF IF (Nn.EQ.l) THEN ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. Auto Correlation ++++++++++++++++ ++++++++++++++++++++++++++++++++++ XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.1O XTOPRIGHT=0.50 YTOPRIGHT=0.90 ELSE IF (Nn.EQ.2) THEN ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 2. Auto Spectra ++++++++++++++++ ++++++++++++++++++++++++++++++++++ INITGRAPHICS=.FALSE. INITFONTS=.FALSE. INITAXES=.TRUE. WINDOWCOLOR=1 IF(BW.EQ.1.0)WINDOWCOLOR=15 XBOTTOMLEFT=0.50 YBOTTOMLEFT=0.10 XTOPRIGHT=1.00 YTOPRIGHT=O.90 ELSE IF (Nn.EQ.3) THEN ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 3. Average Correlation ++++++++++++++++ ++++++++++++++++++++++++++++++++++ WINDOWCOLOR=O IF(BW.EQ.1.0)WINDOWCOLOR=15 XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.10 XTOPRIGHT=0.50 YTOPRIGHT=0.90 ELSE IF (Nn.EQ.4) THEN ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 4. Average Spectra ++++++++++++++++ ++++++++++++++++++++++++++++++++++ WINDOWCOLOR=0 IF(BW.EQ.1.0)WINDOWCOLOR=15 INITGRAPHICS=.FALSE. XBOTTOMLEFT=O.50 YBOTTOMLEFT=O.10 XTOPRIGHT=1.00 YTOPRIGHT=O.90 ELSE IF ((Nn.EQ.5).OR.(Nn.EQ.6)) THEN ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 5. Top Title on the Screen ++++++++++++++++ ++++++++++++++++++++++++++++++++++ INITGRAPHICS=.FALSE. TITLEONLY=.TRUE. INITAXES=.FALSE. INITFONTS=.TRUE. FONTHEIGHT=28 FONTWIDTH=16 WINDOWCOLOR=7 TITLECOLOR=4 IF (BW.EQ.1.0) THEN WINDOWCOLOR=15 TITLECOLOR=4 END IF XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.90 XTOPRIGHT=1.00 YTOPRIGHT=1.00 ELSE IF ((Nn.EQ.7).OR.(Nn.EQ.8)) THEN 0000 000 00000 138 +++++++++++++++++f++++++++++++++++ ++++++++++++++++ 6. Bottom Title Screen ++++++++++++++++ ++++++++++++++++++++++++++++++++++ INITGRAPHICS=.FALSE. TITLEONLY=.TRUE. INITAXES=.FALSE. INITFONTS=.TRUE. FONTHEIGHT=18 FONTWIDTH=9 WINDOWCOLOR=7 TITLECOLOR=0 IF (BW.EQ.1.0) THEN WINDOWCOLOR=15 TITLECOLOR=4 END IF XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.00 XTOPRIGHT=1.00 YTOPRIGHT=0.10 END IF Call plot to print graphs to the screen CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, 1X,Y,AXIS,NUMBER,PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAX, 2XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,EONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,X3,Y3,NUM3,BW 6,PLOT2METHOD) WAIT FOR USER TO HIT RETURN ( only for testing the program , otherwise let the program run w1th no stops) IF ((Nn.EQ.7).OR.(Nn.EQ.8)) THEN READ * *) I END IF IF (Nn.EQ.8) CALL ENDGRAPHICS () RETURN END 139 MYSTAT. FOR C####################################################################### C####################################################################### SUBROUTINE MYSTAT(X,Nn,XMIN,XMAX,NMIN,NMAX,XMEAN,XVAR,XSTD,XVX, 1XSKE,XKUR,XSE,XABS) DIMENSION X(Nn) LOGICAL*1 XABS =-1El4 XMIN=1e14 NMAX=O NMIN= 0 XSUM=0.0 XSUM2=0.0 XSUM3=0.0 XSUM4=0.0 DO 10 I=1,Nn if(xabs) then DO 20 II=1,Nn X(I)=ABS(X(I)) 20 CONTINUE END IF DETERMINE THE MEAN 000 XSUM:XSUM+X(I) XSUM2=XSUM2+X(I)**2.0 FINDING MAXIMUM.AND MINIMUM IE (X(I).GT.XMAX) THEN XMAX=X(I) NMAX=I END IE IE (X(I).LT.XMIN) THEN XMIN=X(I) NMIN=I END IE 10 CONTINUE XN=FLOAT(Nn) XMEAN =XSUM/XN XVAR =XSUM2/XN-XMEAN**2.0 XSTD=SQRT XVAR) XVX=(XSTD XMEAN)*100.0 XSE=XSTD/SQRT(XN) COEFFICIENT OF SKEWNESS AND COEFFICIENT OF KURTOSIS 000 000 DO 30 I=1,Nn XSUM3=XSUM3+(X(I)-XMEAN)**3.0 XSUM4=XSUM4+(X(I)-XMEAN)**4.o 30 CONTINUE XSKE=XSUM3/XSTD**3.0 XKUR=XSUM4/XSTD**4.0 RETURN END $STORAGE:2 $LARGE 140 PLOT.FOR C################################################ C START OF MODULE PLOT FOR ####################### INCLUDE 'EGRAPH.EI' SUBROUTINE PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES, 1INITEONTS,XDATA,YDATA,Y2DATA,NUMBER,PLOTMETHOD,INVERTX,INVERTY, 2XMIN,XMAX,YMIN,YMAX,XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT, 3XTOPRIGHT,YTOPRIGHT,PLOTZCOLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR, 4LABELCOLOR,PLOTCOLOR,TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE, 5TITLEONLY,FONTHEIGHT,FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR, 6X3,Y3,NUM3,BW,PLOT2METHOD) PURPOSE PROGRAMMER REVISION NOTES CALLS GRAPHICSMODE() ON ENTRY LOGICAL*1 LOGICAL*1 LOGICAL*1 LOGICAL*1 LOGICAL*1 REAL*4 REAL*4 INTEGER*2 LOGICAL*1 LOGICAL*1 LOGICAL*1 REAL*4 REAL*4 REAL*4 REAL*4 REAL*4 REAL*4 000000000000000000000000000000000000000000000000000000000000000 TO PLOT A SET OF X-Y DATA USING THE MICROSOFT FORTRAN 5.0 GRAPHICS LIBRARY. PLOT WILL PLOT AN 'UNLIMITED' SET OF X AND Y DATA USING CGA, EGA, VGA OR HERCULES GRAPHICS WITH USER DEFINED PLOT POSITIONING, AXIS LABELLING AND TITLING. LIAM E. GUMLEY, CURTIN UNIVERSITY OF TECHNOLOGY 20 APRIL 1992 BY BASSEM K. KHAFAGI HERCULES GRAPHICS USERS MUST RUN MSHERC.COM (MICROSOFT HERCULES RESIDENT VIDEO SUPPORT ROUTINES, SUPPLIED WITH MICROSOFT FORTRAN 5.0) BEFORE ATTEMPTING TO CALL PLOT. PLOT ALSO USES THE MICROSOFT FORTRAN 5.0 FONT FILE HELVB.FON TO GENERATE CHARACTERS. HELVB.FON MUST BE IN THE CURRENT DIRECTORY WHEN PLOT IS CALLED OTHERWISE THE PROGRAM WILL STOP WITH A 'FONT FILE NOT FOUND' MESSAGE. NOTE THAT THE FONT FILE HELVB.FON ONLY NEEDS TO BE INITIALISED ONCE, ON THE FIRST CALL TO PLOT. THE FONTS REMAIN RESIDENT IN MEMORY FOR SUBSEQUENT CALLS TO PLOT. ALSO NOTE THAT SELECTION OF COLORS MAY NEED TO BE CHANGED ON CGA OR MONOCHROME MONITORS, AS THESE ONLY SUPPORT TWO COLORS (BLACK OR WHITE). THE X AND Y AXIS LABELS ARE WRITTEN IN SCIENTIFIC NOTATION. THE BASE IS WRITTEN NEXT TO THE TICS, AND THE EXPONENT PORTION, TAKEN FROM THE AXIS MAXIMUM VALUE, IS WRITTEN AT THE END OF THE AXIS. SET HIGHEST RESOLUTION ON CGA/EGA/VGA/HGC ALSO CALLS ROUTINES EROM MICROSOFT FORTRAN 5.0 GRAPHICS LIBRARY GRAPHICS.LIB. MAIN PROGRAMS WHICH CALL THIS SUBROUTINE (PLOT) MUST BE LINKED WITH GRAPHICS.LIB. E.G. >FL MAIN.FOR PLOT.FOR GRAPHICS.LIB TRUE : ENTER GRAPHICS MODE FALSE: NO ACTION TRUE : DRAW WINDOW BORDER FALSE: NO ACTION TRUE : ERASE INSIDE WINDOW FALSE: NO ACTION INIT GRAPHICS INIT BORDER INIT WINDOW INIT AXES TRUE : DRAW AXES AND LABELS FALSE: NO ACTION INIT FONTS TRUE : INITIALISE FONT FILE FALSE: NO ACTION XDATA ARRAY OF ORDINATES (ASCENDING) YDATA ARRAY OF COORDINATES NUMBER # OF ELEMENTS IN XDATA, YDATA PLOT METHOD = TRUE : JOIN POINTS WITH LINES = FALSE: DRAW DOTS ATTEBINTS T X = TRUE : X AXIS INVER INVER = FALSE: NO ACTIIHVERTED Y = TRUE : Y AXIS INVERT = FALSE: NO ACTION XMIN MINIMUM TO PLOT FOR X AXIS XMAX MAXIMUM TO PLOT FOR X AXIS YMIN MINIMUM TO PLOT FOR Y AXIS YMAX MAXIMUM TO PLOT FOR Y AXIS X TIC SPACE TIC-LABEL INTERVAL FOR X AXIS Y TIC SPACE TIC-LABEL INTERVAL FOR Y AXIS 0000000000000000000000 10 20 0 00000 00 000 141 REAL*4 X BOTTOM LEFT SCREEN PLOT LOCATION (0.0-1.0) (0.0,0.0) = SCREEN BOTTOM LEFT (l.0,1.0) = SCREEN TOP RIGHT REAL*4 Y BOTTOM LEFT SCREEN PLOT LOCATION (0.0-1.0) REAL*4 X TOP RIGHT SCREEN PLOT LOCATION (0.0-1.0) REAL*4 Y TOP RIGHT SCREEN PLOT LOCATION (0.0-l.0) INTEGER*2 BORDER COLOR COLOR FOR WINDOW BORDER INTEGER*2 WINDOW COLOR COLOR FOR WINDOW BACKGROUND INTEGER*2 AXIS COLOR COLOR FOR PLOT AXES AND TICS INTEGER*2 LABEL COLOR COLOR FOR AXIS LABELS INTEGER*2 PLOT COLOR COLOR FOR PLOTTED DATA INTEGER*2 PLOT2 COLOR COLOR FOR PLOTTED DATA (2ND SET) INTEGER*2 TITLE COLOR COLOR FOR TITLE TEXT CHARACTER*80 X AXIS TITLE TITLE TEXT FOR X AXIS CHARACTER*80 Y AXIS TITLE TITLE TEXT FOR Y AXIS CHARACTER*80 TITLE TITLE TEXT FOR PLOT LOGICAL*1 TITLE ONLY = TRUE : BORDER AND TITLE ONLY = FALSE: NO ACTION INTEGER*2 FONT HEIGHT CHARACTER FONT HEIGHT (PIXELS) INTEGER*2 FONT WIDTH CHARACTER FONT WIDTH (PIXELS) ON EXIT NOTHING CHANGED, NOTHING RETURNED INCLUDE 'FGRAPH.FD' LOGICAL*1 INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lPLOTMETHOD,INVERTX,INVERTY,TITLEONLY,SCRONLY,Y2PLOT,Y3PLOT, 2PLOT2METHOD INTEGER*2 DUMMY,XWIDTH,YHEIGHT,IXl,IY1,IX2,IY2,BORDERCOLOR, 1WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR,TITLECOLOR,NUMBER,I,IX, 2IY,TICLENGTH,PLOTZCOLOR,FONTHEIGHT,FONTWIDTH,TITLELENGTH, 3PLOT3COLOR REAL*4 X,Y,XMIN,YMIN,XMAX,YMAX,XDATA(NUMBER),YDATA(NUMBER), lY2DATA(NUMBER),XWl,YW1,XW2,YW2,XTICSPACE,YTICSPACE,XBOTTOMLEFT, ZYBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT,XT,YT,X3(NUM3),Y3(NUM3) CHARACTER TICLABEL*8,FONTFILE*12,FONTSTRING*7,TITLE*80, lXAXISTITLE*80,YAXISTITLE*80,FONTCOMMAND*10,TITLE1*80,TITLE2*80, 2TITLE3*80,TITLE4*80,TITLE5*80,TITLE6*80,TITLE7*80,TITLE8*80 RECORD / VIDEOCONFIG / SCREEN RECORD / WXYCOORD / WXY RECORD / XYCOORD / POSITION COMMON SCREEN SET THE NAME OE THE EONT FILE To USE (MUST BE IN CURRENT DIR) FONTFILE='HELVB.FON' FONTCOMMAND="T'HELV'" SET THE GRAPHICS MODE To HIGHEST POSSIBLE RESOLUTION IE (INITGRAPHICS) CALL GRAPHICSMODE () REGISTER FONT FOR LABELS IE (INITFONTS) THEN CALL UNREGISTEREONTS () IE (REGISTEREONTS(FONTFILE).LT.0) THEN WRITE (*,10) FONTFILE EORMAT (' EONT FILE ',A12,' NOT FOUND IN PLOT') STOP END IE SET UP EONT FOR USE WRITE (FONTSTRING,20) FONTHEIGHT,FONTWIDTH EORMAT ('H',I2.2,'W',IZ.2,'B') DUMMY=SETEONT(FONTCOMMAND//FONTSTRING) END IE GET SCREEN DIMENSIONS XWIDTH=SCREEN.NUMXP;§§EES YHEIGHT=SCREEN.NUMY SET A.VIEW PORT, WITH CORNER COORDINATES DEFINED AS A FRACTION OF THE TOTAL SCREEN SIZE (0.0,0.0) IS BOTTOM LEFgFogcgggfiEN . 1.0 IS TOP RIGHT TRANSLATE CORNER COORDINATES TO PIXEL COORDINATES Ix1=INT(XBOTTOMLEET*FLOAT(XWIDTH)) IX2=INT(XTOPRIGHT*FLOAT(XWIDTH)) IYl=YHEIGHT—INT(YTOPRIGHT*FLOAT(YHEIGHT)) IY2=YHEIGHT-INT(YBOTTOMLEFT;§L%$§(YHEIGHT)) TVIEWPORT 1x1 IY , , CACREATE.A REAL COORDINATE WINDOW, WITH GRAPH AREA SIZED AT 60% OF THE TOTAL VIEWPORT SIZE (LEAVES 15% EITHER SIDE FOR LABELS AND TICS) SCALE=0.60 IF THE GRAPH OCCUPIES LESS THAN 50% OF THE SCREEN IN EITHER THE X OR Y DIRECTION, THEN THE LABELS WON'T FIT. 00 000 000 000 0000 0000 30 40 142 80 IF THE GRAPH IS LESS THAN 50% OF SCREEN IN X OR Y, THEN SIZE THE GRAPH AREA AT 30% OF THE VIEWPORT SIZE. IF (XTOPRIGHT-XBOTTOMLEFT.LT.0.5.0R.YTOPRIGHT-YBOTTOMLEFT.LT.0.5) ISCALE=0.30 IF (SCRONLY) THEN XW1=XMIN YW1=YMIN XW2=XMAX YW2=YMAX GO TO 30 END IF XW1=XMIN-0.5*(1.0-SCALE)*ABS(XMAX-XMIN) YWl=YMIN-0.5*(l.0-SCALE)*ABS(YMAX-YMIN) XW2=XMAX+0.5*(1.0-SCALE)*ABS(XMAX-XMIN) YW2=YMAX+0.5*(1.0-SCALE)*ABS(YMAX-YMIN) DUMMY=SETWINDOW(.TRUE.,XW1,YW1,XW2,YW2) DRAW A BORDER AROUND THE VIEWPORT IF (INITBORDER) THEN IF (INITWINDOW) THEN DUMMY=SETCOLOR(WINDOWCOLOR) DUMMY=RECTANGLE WK$GFILLINTERIOR,XW1,YW1,XW2,YWZ) END IF _ DUMMY=SETCOLOR(BORDERCOLOR) DUMMY=RECTANGLE_W($GBORDER,XWl,YW1,XW2,YWZ) END IF SKIP TO TITLE DRAWING IF REQUIRED IF (TITLEONLY) GO TO 90 IF (SCRONLY) GO TO 40 IF (.NOT.SCRONLY) GO TO 50 This is the screen design of the openning logo TITLELENGTH=80 DUMMY=SETCOLOR(8) IF (BW.EQ.1.0)DUMMY=SETCOLOR(15) YYW1=YW1+0.42 YYW2=YW2-0.42 DUMMY=RECTANGLE W($GFILLINTERIOR,XW1,YYW1,XW2,YYWZ) DUMMY=SETCOLOR(I2) IF (BW.EQ.1.0)DUMMY=SETCOLOR(0) DUMMY=RECTANGLE;W($GBORDER,XWl,YYW1,XW2,YYW2) POSITION THE TITLE CENTERED (MORE OR LESS) ABOVE GRAPH YT=0.5*(YMAX+YMIN) CALL MOVETO W (1.0*(XMAX+XMIN),YT WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-(TITLELENGTH*FONTWIDTH)/2 IY=IY—FONTHEIGHT/2 enter the program name (Title: from the main program) DUMMY=SETCOLOR(0) CALL MOVETO (IX,IY,POSITION) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) DUMMY=SETCOLOR(11) IF (BW.EQ.1.0)DUMMY=SETCOLOR(0) CALL MOVETO (IX,IY+2,POSITION) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) Standard logo TITLE1=' DEPARTMENT OF CIVIL ENGINEERING' TITLE2=: ' MICHIGAN STATE UNIVERSITY' TITLE3= TITLE4=' LOAD RESEARCH' TITLE5=' PROGRAMS FOR DISSERTATION' TITLE6=' ' TITLE7=' BASSEM K. KHAFAGI' TITLE8=' OCTOBER 1992' Drawing the screen 000 50 0000 C60 61 70 143 DUMMY=SETCOLOR(15) CALL MOVETO (IX,IY-200,POSITION) CALL OUTGTEXT (TITLEl(l:TITLELENGTH)) CALL MOVETO (IX,IY-160,POSITION) CALL OUTGTEXT (TITLE2(1:TITLELENGTH)) CALL MOVETO (IX,IY-120,POSITION) CALL OUTGTEXT (TITLE3(1:TITLELENGTH)) CALL MOVETO (IX,IY-80,POSITION) CALL OUTGTEXT (TITLE4(1:TITLELENGTH)) CALL MOVETO (IX,IY+80,POSITION) CALL OUTGTEXT (TITLE5(1:TITLELENGTH)) CALL MOVETO (IX,IY+120,POSITION) CALL OUTGTEXT (TITLE6(1:TITLELENGTH)) CALL MOVETO (IX,IY+160,POSITION) CALL OUTGTEXT (TITLE7(1:TITLELENGTH)) CALL MOVETO (IX,IY+200,POSITION) CALL OUTGTEXT (TITLE8(1:TITLELENGTH)) Adding shadow to the text DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY-198,POSITION) CALL OUTGTEXT (TITLE1(1:TITLELENGTH)) DUMMY=SETCOLOR(1) CALL MOVETO (IX,IY-158,POSITION) CALL OUTGTEXT (TITLE2(1:TITLELENGTH)) DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY-118,POSITION) CALL OUTGTEXT (TITLE3(1:TITLELENGTH)) CALL MOVETO (IX,IY—79,POSITION) CALL OUTGTEXT (TITLE4(1:TITLELENGTH)) CALL MOVETO (IX,IY+84,POSITION) CALL OUTGTEXT (TITLE5(1:TITLELENGTH)) CALL MOVETO (IX,IY+122,POSITION) CALL OUTGTEXT (TITLE6(1:TITLELENGTH)) DUMMY=SETCOLOR(1) CALL MOVETO (IX,IY+162,POSITION) CALL OUTGTEXT (TITLE7(1:TITLELENGTH)) DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY+202,POSITION) CALL OUTGTEXT (TITLE8(1:TITLELENGTH)) SKIP THE AXES DRAWING IF NOT NEEDED IF (.NOT.INITAXES) GO To 140 DRAW AXES DUMMY=SETCOLOR(AXISCOLOR) CALL MOVETO W (XMIN,YMAX,WXY) DUMMY=LINETO‘W(XMIN,YMIN) DUMMY=LINETO‘W(XMAX,YMIN) DUMMY=LINETO‘W(XMAX,YMAX) DUMMY=LINETO_W(XMIN,YMAX) PUT TIC MARKS AND LABELS ON AXES AT USER DEFINED SPACES TIC MARKS ARE ALWAYS 2 PIXELS LONG TICLENGTH=2 TIC MARKS AND LABELS EOR X AXIS DO 70 X=XMIN,XMAX,XTICSPACE INVERT X AXIS PLOTTING IE REQUIRED XT=X IE (INVERTX) XT=XMAX-ABS(XT—XMIN) CALL MOVETO W (XT,YMIN,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD DUMMY=SETCOLOR(AXISCOLOR) DUMMY=LINETO(IX,IY+TICLENGTH) TIC LABEL IS 5 CHARACTERS LONG WRITE (TICLABEL,60) X EORMAT (1PE8.1) EORMAT (1F5.1) EORMAT (1F6.2) MAKE ROOM FOR THE TIC LABEL IX=IX—(5*FONTWIDTH)/2 IY=IY+(FONTHEIGHT/2) CALL MOVETO (IX,IY,POSITION) DUMMY=SETCOLOR(LABELCOLOR) CALL OUTGTEXT (TICLABEL(1:5)) CONTINUE PUT TIC MARK AT XMAX IN CASE IT WAS MISSED 80 90 100 110 144 XT=XMAX IE (INVERTX) XT=XMIN CALL MOVETO W (XT,YMIN,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD DUMMY=SETCOLOR(AXISCOLOR) DUMMY=LINETO(IX,IY+TICLENGTH) TIC MARKS AND LABELS FOR Y AXIS DO 80 Y=YMIN,YMAX,YTICSPACE YT Y INVERT Y AXIS PLOTTING IE REQUIRED IF (INVERTY) YT=YMAX-ABS(YT-YMIN) CALL MOVETO W (XMIN,YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD DUMMY=SETCOLOR(AXISCOLOR) DUMMY=LINETO(IX—TICLENGTH,IY) TIC LABEL IS 5 CHARACTERS LONG WRITE (TICLABEL,61) Y MAKE ROOM FOR THE TIC LABEL IX=IX-(6*FONTWIDTH) IY=IY-(FONTHEIGHT/2) CALL MOVETO (IX,IY,POSITION) DUMMY=SETCOLOR(LABELCOLOR) CALL OUTGTEXT (TICLABEL(1:6)) CONTINUE PUT TIC MARK AT YMAX IN CASE IT WAS MISSED YT=YMAX IF (INVERTY) YT=YMIN CALL MOVETO W (XMIN,YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD DUMMY=SETCOLOR(AXISCOLOR) DUMMY=LINETO(IX-TICLENGTH,IY) WRITE THE GRAPH TITLE AT TOP OF GRAPH REGION DUMMY=SETCOLOR(TITLECOLOR) POSITION THE TITLE CENTERED (MORE OR LESS) ABOVE GRAPH TITLELENGTH=80 IF (TITLE(TITLELENGTH:TITLELENGTH).EQ.' ') THEN TITLELENGTH=TITLELENGTH-1 GO TO 100 END IF YT=YMAX IE (TITLEONLY) YT=0.5*(YMAX+YMIN) CALL MOVETO W (0.5*(XMAX+XMIN),YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-(TITLELENGTH*FONTWIDTH)/2 IF (TITLEONLY) THEN IY=IY—EONTHEIGHT/2 ELSE IY=IY-2*FONTHEIGHT END IF CALL MOVETO (IX,IY,POSITION) DUMMY=SETCOLOR(TITLECOLOR) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) Go To RETURN STATEMENT IF ONLY TITLE TO BE DRAWN IE (TITLEONLY) GO To 170 WRITE THE XOAXIS TITLE TITLELENGTH= IF (XAXISTITLE(TITLELENGTH:TITLELENGTH).EQ.‘ ') THEN TITLELENGTH=TITLELENGTH-1 Go TO 110 END IF POSITION THE TITLE CENTERED (MORE OR LESS) BELOW GRAPH CALL MOVETO W (0.5*(XMAX+XMIN),YMIN,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-(TITLELENGTH*FONTWIDTH)/2 IY=IY+(3*FONTHEIGHT)/2 CALL MOVETO (IX,IY,POSITION) DUMMY=SETCOLOR(TITLECOLOR) 000 120 130 140 150 145 CALL OUTGTEXT (XAXISTITLE(1:TITLELENGTH)) WRITE THE Y AXIS TITLE VERTICALLY, SINCE WE CAN'T ROTATE TITLELENGTH=80 IE (YAXISTITLE(TITLELENGTH:TITLELENGTH).EQ.' ') THEN TITLELENGTH=TITLELENGTH-l GO TO 120 END IF POSITION THE TITLE CENTERED (MORE OR LESS) TO LEFT OF GRAPH CALL MOVETO W (XMIN,0.5*(YMAX+YMIN),WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-8*FONTWIDTH IY=IY-(TITLELENGTH*FONTHEIGHT)/2 DUMMY=SETCOLOR(TITLECOLOR) NOW WRITE THE CHARACTERS ONE AT A.TIME VERTICALLY Do 130 I=1,TITLELENGTH CALL MOVETO (IX,IY,POSITION) CALL OUTGTEXT (YAXISTITLE(I:I)) IY=IY+EONTHEIGHT CONTINUE PLOT THE DATA.IN XDATA AND YDATA DUMMY=SETCOLOR(PLOTCOLOR) DEFINE CLIPPING REGION GET TOP LEFT POSITION CALL MOVETO W (XMIN,YMAX,WXY) CALL GETCURRENTPOSITION (POSITION) CALL GETPHYSCOORD (POSITION.XCOORD,POSITION.YCOORD,POSITION) IX1=POSITION.XCOORD IY1=POSITION.YCOORD GET BOTTOM RIGHT POSITION CALL MOVETO W (XMAX,YMIN,WXY) CALL GETCURRENTPOSITION (POSITION) CALL GETPHYSCOORD (POSITION.XCOORD,POSITION.YCOORD,POSITION) IX2=POSITION.XCOORD IY2=POSITION.YCOORD SET UP NEW VIEWPORT FOR CLIPPING CALL SETVIEWPORT (IX1,IY1,IX2,IY2) DUMMY=SETWINDOW(.TRUE.,XMIN,YMAX,XMAX,YMIN) PLOT THE FIRST POINT, INVERTING IE NECESSARY XT=XDATA(1) YT=YDATA(1) IE (INVERTX) XT=XW2-ABS(XT—XW1) IE (INVERTY) YT=YW2-ABS(YT—XW2) CALL MOVETO W (XT,YT,WXY) IE (.NOT.PLOTMETHOD) DUMMY=SETPIXEL W(XT,YT) PLOT THE REST OF THE DATA, INVERTING IF NECESSARY Do 150 I=2,NUMBER XT=XDATA(I) YT=YDATA(I) IF (INVERTX) XT=XW2—ABS(XT-XW1) IF (INVERTY) YT=YW2—ABS(YT—YW1) IF (PLOTMETHOD) THEN DUMMY=LINETO_W(XT,YT) ELSE CALL MOVETO W (XT,YT,WXY) DUMMY=SETPIXEL;W(XT,YT) END IF CONTINUE PLOT THE SECOND SET OF DATA ON THE SAME X-Y AXIS IE(.NOT.Y2PLOT) Go TO 170 DUMMY=SETCOLOR(PLOTZCOLOR) XT=XDATA(1) YT=Y2DATA(1) IF (INVERTX) XT=XW2-ABS(XT—XW1) IE (INVERTY) YT=YWEEAaiéYT-XW2) CALL MOVETO W XT, , IF (.NOT.PLOT2METHOD) DUMMY=SETPIXEL W(XT,YT) PLOT THE REST OF THE DATA, INVERTING IE NECESSARY DO 160 I=2,NUMBER XT=XDATA(I) YT=Y2DATA(I) IF (INVERTX) XT=XW2-ABS(XT-XW1) IE (INVERTY) YT=YW2-ABS(YT-YW1) IF (PLOT2METHOD) THEN 160 000 165 C c 170 146 DUMMY=LINETO W(XT,YT) ELSE - CALL MOVETO W (XT,YT,WXY) DUMMY=SETPIXEL;W(XT,YT) END IF CONTINUE PLOT THE THIRD SET OF DATA ON THE SAME X-Y AXIS IE(.NOT.YSPLOT) GO TO 170 DUMMY=SETCOLOR(PLOT3COLOR) XT=X3(1) YT=Y3(I) IE (INVERTX) XT=XW2~ABS(XT-Xw1) IF (INVERTY) YT=YW2-ABS (YT-XW2) CALL MOVETO W (XT,YT,WXY) IE (.NOT.PLOT2METHOD) DUMMY=SETPIXEL W(XT,YT) PLOT THE REST OF THE DATA, INVERTING IF NECESSARY DO 165 I=2,NUM3 XT=X3(1) YT=Y3(I) IF (INVERTX) XT=Xw2-ABS(XT—Xw1) IF (INVERTY) YT=YW2—ABS(YT-YW1) IE (PLOT2METHOD) THEN DUMMY=LINETO W(XT,YT) ELSE ‘ CALL MOVETO W (XT,YT,WXY) DUMMY=SETPIXEL W(XT,YT) END IF — CONTINUE ALL DONE CALL UNREGISTERFONTS () RETURN END C#############################??######################################## 000 000 SUBROUTINE GRAPHICSMODE INCLUDE 'FGRAPH.FD' INTEGER*2 DUMMY,MAXX,MAXY RECORD /VIDEOCONFIG/ MYSCREEN COMMON MAXX,MAXY FIND GRAPHICS MODE. CALL GETVIDEOCONFIG (MYSCREEN) SELECT CASE( MYSCREEN.ADAPTER ) CASE( $CGA ) DUMMY=SETVIDEOMODE($HRESBW) CASE( $OCGA ) DUMMY=SETVIDEOMODE($ORESCOLOR) CASE( SEGA, $OEGA ) IF (MYSCREEN.MONITOR.EQ.$MONO) THEN DUMMY=SETVIDEOMODE($ERESNOCOLOR) ELSE DUMMY=SETVIDEOMODE($ERESCOLOR) END IF CASE( $VGA, $OVGA, $MCGA ) IE (MYSCREEN.MONITOR.EQ.$MONO) THEN DUMMY=SETVIDEOMODE($VRESZCOLOR) ELSE DUMMY=SETVIDEOMODE($VRE316COLOR) END IF CASE( $HGC ) DUMMY=SETVIDEOMODE($HERCMONO) CASE DEFAULT DUMMY=0 END SELECT IF (DUMMY.EQ.0) STOP 'ERROR: CANNOT SET GRAPHICS MODE' DETERMINE THE MINIMUM.AND MAXIMUM DIMENSIONS. CALL GETVIDEOCONFIG (MYSCREEN) MAXX=MYSCREEN.NUMMPIXELS-l MAXY=MYSCREEN.NUMYPIXELS-l END C############################?f######################################### SUBROUTINE ENDGRAPHICS INCLUDE 'FGRAPH.FD' 147 INTEGER*2 DUMMY DUMMY=SETVIDEOMODE($DEFAULTMODE) RETURN END *1-*I-**)f*fi-**if*i-**’+**'**>f*i-**'**it*fl-**Jf*i-**16*fl'*1-** 0000000*************#**#*********t******#*********** 0000000000 148 FLOOR.FOR * * * k * * * * * * i * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * k * * * * * * * * FLOOR.FOR * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * THIS PROGRAM READS THE VOLTAGE DATA COLLECTED FROM.THE PLATE FORM AND FINDS THE CORRESPONDING FORCES. THE DATA IS STORED AT A.RATE OF 100Hz BY 'TRAP.BAS'. THIS PROGRAM USED WITH ANY PERIODIC JUMPING TEST. HOWEVER, THE DEFAULT IS NOW 2 HZ PERIODIC JUMPING. IF THE FREQUENCY DIFFERS FROM 2 HZ, YOU MAY NEED TO CHANGE THE COUNTER IN DO LOOP 80. THE VALUE OF THE COUNTER IN THIS LOOP (30) DETEMINES THE RANGE OF SEARngNG FOR THE FIRST MAJOR PEAK. SEE THE COMMENT JUST BEFORE LOOP . THE PROGRAM NOW DOES NOT INCLUDE THE MASS OF THE PARTICIPENTS IN THE CALCULATION OF THE MASS, DAMP, AND TOTAL FORCE MATRICES. IF YOU NEED TO CONSIDER WITH THE MASS OF THE FLOOR PLUS PART., YOU MAY HAVE TO CHANGE THE MATRIX FMASS(I) IN LOOP 250 TO M(I). ALSO, YOU WILL NEED TO CHANGE THE DEFINITON OF MASS MA IN LOOP 360 TO BE EQUAL TO MASS(I,J) INSTEAD OD TMASS(I,J) #* #* ** ** ** ** +* #* ** ** ** ** ** ** ** ** ** ** *1 ** *w #* ** ** ** ** #fi ** ** ** ** t¥ #* OCTOBER 4, 1989 BASSEM K. KHAFAGI YING X. CAI UNIVERSITY OF IDAHO MOSCOW, ID 83843 *fi'**lf*I-*X-*fi'**=f**'**J&*X-**>+*I-*I-**l$*fl'**!f*I-**I¥*1-** **'**'**'*i-*#'**-*i-*I-*i-*#~*fl-**-*fl-ii-*I~*I-¥*'** *I-#i-*1-**-* U.S.A * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * k * * * * * * * * * * * * * * * * * * * * * * * * * * * * i * * * * * * ***************************************************************** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA.INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ PROGRAM FLOOR DIMENSION VOLT(1000,9),X(300,9),CALI(9),D(150,9),DEF(128,9) DIMENSION A(128,9),B(128,9),VEL(128,9),ACCL(128,9),EIT(128,9) DIMENSION PHIT(9),MN(9),MNN(9),MA(9,9),DA(9,9),DAA(9,9),DAM(9,9) DIMENSION DAMP(9,9),VNODE(9),ANODE(9),XNODE(9),FA(9),FV(9),FX(9) DIMENSION STIF(9,9),EORCES(128,9),PART(40),FMASS(9),TOTEORCE(124) DIMENSION TMASS(9,9),BVOLT(1,9) I REAL*8 K(9,9),M(9) LAMBDA(9),E(9),PHI(9,9),MASS(9,9),NODE(9) REAL MN,MNN,MA CHARACTER FNAME*12,DATE*12,STIME*12,ETIME*12,FORCE*12,TFORCE*12 CHARACTER FNAMEl*3,FNAME2*1,FNAME3*1,FNAME4*1,FNAME5*6,REPORT*12 LOGICAL WANTX NOW, WE READ THE FILE NAME OF THE DATA SET. WE WILL USE THIS FILE NAME AS A BASE FOR NAMING THE OUTPUT FILES. IF THE DATA FILE NAME IS 20PJ2-10.0C9, IT MEANS THAT THIS DATA SET FOR 20 PEOPLE PERFORMING PERIODIC JUMPING (PJ) AND THIS IS GROUP (2) AND TEST NO. (10) THAT WAS DONE IN OCT. 1989 (0C9). THE OUTPUT FILES WILL BE IN THE FORM OF 20PF2F10.0C9, 20PF2T10.0C9, AND 20PF2R10.0C9. THE PART (PFZF,PF2T,PF2R) REPRESENT PERIODIC JUMPING (P), FINAL FORCES (F) FOR GROUP (2) THEN THE LAST LETTER IN THIS PART REPRESENT THE FILE CONTENT. (F) REPRESENTS NODAL FORCES FOR THIS TEST, (T) REP- 000 0 000000000000 000000 0000000000 Tocoqmmawm 12 10 15 18 19 210 211 212 149 RESENTS TOTAL FORCES FROM THE PERIODIC JUMPING ON THE FLOOR, AND (R) REPRESENTS THE FINAL REPORT FOR THIS TEST. PRINT 2 READ(*,3)NODEBASE PRINT 5 READ(*,4)FNAME1,FNAME2,ENAME3,FNAME4,ENAME5 WRITE(ENAME,6)ENAME1,ENAMEz,ENAME3,FNAME4,FNAME5 WRITE(EORCE,7)ENAME1,ENAME3,FNAME5 WRITE(TEORCE,8)ENAME1,FNAME3,FNAME5 WRITE(REPORT,9)ENAME1,FNAME3,FNAME5 EORMAT(1X,' WHICH NODE WILL BE USED AS A.BASE FOR CALCULATIONS') EORMAT(Il) FORMAT(A3,A1,A1,A1,A6) FORMAT(1X,' WHAT 13 NAME OF INPUT FILE?') EORMAT(A3,A1,A1,A1,A6) FORMAT(A3,'F',A1,'F',A6) FORMAT(A3,'E',A1,'T',A6) FORMAT(A3,'F',A1,'R',A6) OPEN(1,FILE=ENAME) OPEN(2,EILE='INPUT') OPEN(3,FILE=EORCE) OPEN(4,FILE=TEORCE) OPEN(5,FILE=REPORT) OPEN(6,FILE='PHI') +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 2. REDUCING DATA TO 128 SAMPLE +++++++++++++ +++++++++++++++++++++++++++++++++++++++ READ THE DATE, STARTING TIME, AND ENDING TIME OF THE TEST. THEN, DELETE THE FIRST FOUR-COLUMNS OF DATA.(STRAIN GAGE DATA), AND SAVE THE LVDT VOLTAGE IN THE FILE "LVDT.VLT" READ(l,12)DATE,STIME,ETIME EORMAT(A12//,A12/,A12) DO 10 I=1,1000 READ(1,15) Y,Y,Y,Y, (VOLT(I,J) ,J=1, 9) CONTINUE FORMAT(13F8.0) SELECT THE FIRST VOLTAGE READING AS THE BASE TO NORMALIZE THE VOLTAGE READING AGAINST THE PART. WEIGHTS AND SAVE IT IN MATRIX "BVOLT" DO 18 I=1,9 BVOLT(1,I)=VOLT(1,I) CONTINUE SEARCH THE FIRST MAJOR VOLTAGE PEAK VALUE AND SELECT 300 SAMPLE POINTS SUCH THAT 149 READINGS BEFORE IT AND 150 AFTER. THE SEARCH WILL USE LVDT NO 5 AS A BASIS SINCE IT IS A.MID. NODE IN THE FLOOR SYSTEM. .ALSO IF NVMAX IS >850 WE WILL ASSUME IT 850 AND IF IT IS < 150 WE WILL ASSUME IT 150. IF THE TEST ENDED BEFORE 1000 SAMPLES THE SEARCH WILL BE TERMINATED AT "NEND" VARIABLE J=NODEBASE N=0 NEND=1000 DO 19 I=1,1000 N=N+1 IF((VOLT(I,J).EQ.0.0).AND.(VOLT(I+1,J+1).EQ.0.0)) GO TO 210 CONTINUE NEND=N WRITE(*,212) NEND FORMAT(1X,'THE TEST ENDED AT SAMPLE NO: (',I4,')‘) VMAX=0.0 N=0 DO 20 I=1,1000 N=N+1 0000 0000000000 000000000000 00000 20 LAMA 80 150 IF (ABS (VOLT (I, J) ) .GT.VMAX) THEN VMAX=ABS(VOLT(I,J)) NVMAX=N Go To 20 ENDIF CONTINUE =NVMAX MEND=NEND—150 IF(NVMAX.GE.MEND)THEN NVMAX=MEND ELSE IF(NVMAX.LE.150)THEN NVMAX=150 END IF DO 22 I=1,300 DO 25 J=1,9 VOLT(I,J)=VOLT(NVMAX+I-150,J) CONTINUE CONTINUE NORMALIZE DATA BY SUBTRACTING EACH RAW OF DATA.FROM THE FIRST VLOTAGE READING WHICH REPRESENT THE STATIC WEIGHT OF THE FLOOR AND PARTICIPANTS. THEN THE DEFLECTION IS CALCULATED BY DIVIDING EACH COLUMN BY THE CORESPONDING CALIBRATION FACTOR OF THAT LVDT. EACH COLUMN REPRESENT THE DATA READ BY ONE OF THE LVDT'S. THE CALIBRATION FACTORS ARE STORED IN THE FIRST PART OF FILE "INPUT" DO 30 J=1,9 READ(2,*) CALI(J) DO 40 I=1,300 X(I,J)=(VOLT(I,J)—BVOLT(1,J))/(-3276.7*CALI(J)) CONTINUE CONTINUE THE DEFLECTION WILL BE REDUCED NOW TO 150-SAMPLE BY DELETING EVERY OTHER SAMPLE, AND SAVE AS MATRIX D(150,9) DO 60 I=1,300,2 N=(I+1)/2 DO 70 J=1,9 D(N,J)=X(I,J) CONTINUE CONTINUE SEARCH THE FIRST MAJOR DEFLECTION PEAK VALUE AND SELECT 128 SAMPLE POINTS SUCH THAT 4 READINGS BEFORE IT AND 123 AFTER. THE SEARCH WILL USE LVDT NO 5 AS A BASIS SINCE IT IS A MID. NODE IN THE FLOOR SYSTEM. SINCE THE FREQ. OF THE TEST IS 2 HZ, AND THE RATE OF COLLECTING DATA IS NOW REDUCED TO 50 HZ, A FULL JUMP WILL BE REPRESENTED BY 25 SAMPLES. IN ORDER TO BE ABLE TO GET TO THE MAX. PEAK, WE MUST SEARCH FOR IT IN A.RANGE GREATER OR EQUAL THAN 25. SO, WE SET THE LOOP COUNTER TO 30. J=NODEBASE DMAX=0.0 NMAX=0 N=o DO 80 I=1,3o N=N+1 IE (D(I,J).GT.DMAX) THEN DMAX=D(I,J) NMAX=N GO To 80 ENDIF CONTINUE MAX=NMAX NOW THE REDUCED SAMPLE OF DEFLECTIONS WILL BE STORED IN MATRIX DEF(128,9) TO USE IT IN THE PROGRAM. THE REDUCED NO. OF SAMPLES (128) WILL BE DENOTED AS "NS" NS=128 DO 90 I=1,NS DO 100 J=1,9 0000000000000 000 000000000000 0000000000 151 DEF(I,J)=D(NMAX+I-S,J) 100 CONTINUE 9o CONTINUE +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 3. FINDING FITTED DEFLECTION +++++++++++++ +++++++++++++++++++++++++++++++++++++++ PERFORM PERIODIC REGRESSION OF SELECTED DEFLECTION SAMPLES BY USING THE FOURIER SERIES TO GET FITTED DEFLECTION. CALCULATE FOURIER COEFFICIENTS PI=3.141593 DO 110 J=1,9 SUM1=0. DO 120 I=1,NS SUM1=SUM1+DEE(I,J) 120 CONTINUE .A(0,J)=SUM1/NS DO 130 I=1,20 SUM2=O SUM3=O DO 140 IJ=1,NS SUM2=SUM2+DEF(IJ,J)*COS(2*I*PI*IJ/NS) SUM3=SUM3+DEF(IJ,J)*SIN(2*I*PI*IJ/NS) 14o CONTINUE .A(I,J)=2*SUM2/NS B(I,J)=2*SUM3/NS 13o CONTINUE CALCULATE PREDICTED DEFLECTION HISTORY D0 150 I=1,NS SUM1=0 SUM2=0 DO 160 JJ=1,20 SUM1=SUM1+A(JJ,J)*COS(2*PI*JJ*I/NS) SUM2=SUM2+B(JJ,J)*SIN(2*PI*JJ*I/NS) 160 CONTINUE FIT(I,J)=A(0,J)+SUM1+SUM2 150 CONTINUE 110 CONTINUE +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 4. VELOCITY AND ACCELERATION +++++++++++++ +++++++++++++++++++++++++++++++++++++++ ************‘k**************************************************** PROGRAM TO DO CENTRAL DIFFERENTIAL NOTE: H IS THE TIME INTERVAL, H= l/(FREQUENCE OF SELECTED SAMPLE POINTS). SINCE THE SAMPLES HAVE BEEN REDUCED TO 50 HZ, H=.02 SEC. H=0.02 DO 180 J=1,9 DO 190 I=3,126 VEL(I,J)=(-§IT(I+2,J)+8*FIT(I+1,J)-8*FIT(I-l,J)+FIT(I-2,J)) 2 H ( ) ACCL(I,J)=(-FIT(I+2,J)+16*FIT(I+1,J)-30*FIT(I,J)+16*FIT(I-1 ,J)-FIT(I-2,J))/(12*H*H) 190 CONTINUE 180 CONTINUE +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 5. FLOOR AND PART. MASS +++++++++++++ +++++++++++++++++++++++++++++++++++++++ HERE, WE CALCULATE THE MASS MATRIX FOR THE FLOOR SYSTEM AND THE PARTICIPANTS USING THE WEIGHT AND LOCATION OF EACH PARTICIPANT AND APPLYING IT TO THE CORRECT NODE. 00000000000 220 000 * 000 225 000 226 000 227 0000 240 000 152 THE ARRAYS INCLUDE: PART(40) - WEIGHTS AND LOCATIONS OF PARTICIPANTS M(9) - MASSES OF PARTICIPANTS BROKEN UP INTO 9 NODES NODE(9) - PARTICIPANT AND FLOOR MASS FOR 9 NODES MASS(9,9) - MASS MATRIX TO BE SAVED IN FILE "MASS" READ THE PARTICIPANT WEIGHT-LOCATION FILE, (THIS IS IN THE SECOND PART OF THE FILE "INPUT". N=0 DO 220 J=1, 4o READ(2, *) PART(J PART(J) =PART(J)/386. 4 N=N+1 IF(PART(J).GT.0.0) THEN NPART=N GO To 220 ENDIF CONTINUE CALCULATE TOTAL WEIGHT CORRESPONDING TO EACH NODE NODE(1)— =P?RT$32£;§ART(29)+2. ./3. *(PART(37)+PART(27)+PART(19)) T NODE(2)- =1. /3. *(PART(19)+PART(20))+2. /3. *(PART(7)+PART(8)) 2./9. *(PART(17)+PART(18))+PART(9)+PART(10) NODE(3)=PART(30)+PART(40)+2. /3. *(PART(20)+PART(28)+PART(38)) +4. /9. *PART(18) NODE(4) =PART(35)+PART(25)+2. /3. *PART(15)+2. /9. *(PART(17)+PART(13)) +1./3 .*(PART(37)+PART(27)+PART(33)+PART(23)) NODE(S) =1. /9. *(PART(17)+PART(18)+PART(13)+PART(14)) +1 ./3. *(PART(7)+PART(8)+PART(3)+PART(4)+PART(15)+PART(16)) +PART(5)+PART(6) NODE(6) =PART(26)+PART(36)+2. /3. *PART(16)+2. /9. *(PART(18)+PART(14)) +1. /3. *(PART(28)+PART(38)+PART(24)+PART(34)) :NODE(7) =PART(31)+PART(21)+2. /3. *(PART(33)+PART(23)+PART(11)) +4 /9. *PART(13) NODE(8)- =PART(1)+PART(2)+2. /3. *(PART(3)+PART(4)) +2 9.*(PART(13)+PART(14))+1 /3. *(PART(11)+PART(12) NODE(9) =PART(22)+PART(32)+2. /3. *(PART(24)+PART(3 4)+PART(12) +4. /9. *PART(14) ) ) INPUT THE MASS OF THE FLOOR ALONE IN ARRAY FMASS(9) DO 225 I=1,9 FMASS(I)=0.358592 CONTINUE FMASS(2)=0.222567 FMASS(5)=O.222567 FMASS(8)=0.222567 PUT FLOOR MASSE IN THE FORM OF A 9X9 FLOOR MASS MATRIX DO 226 I= 1, 9 J=I CONTINUE ADD MASS OF PARTICIPANTS TO MASS OF FLOOR FOR EACH NODE DO 227 I= 1, 9 M(I)=NODE(I)+FMASS(I) CONTINUE TMASS(I,J)=FMASS(J) PUT NODAL MASSES IN THE FORM OF A 9X9 MASS MATRIX DO 230 JI=1, 9 =I MASS(I,J)=M(J) CONTINUE FORMAT(9F10.6) 00000000000000000000000000000000000 00000 000 000 000 000 000 000 153 +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 6. DAMPING MATRIX +++++++++++++ +++++++++++++++++++++++++++++++++++++++ PROGRAM TO SET UP AND SOLVE CHARACTERISTIC VALUE PROBLEMS OF THE TYPE: (LAMBDA)*M*X = K*X THIS PROGRAM READS THE MATRICES K AND M AND CONVERTS IT TO THE STANDARD EIGENVALUE FORM, (LAMBDA)*X = AFK BY THE FOLLOWING TRANSFORMATION: (L.AMBDA)*(M**0. 5)*X = (M**- -.0 5)*K*(M**- -.0 5)*(M**0. 5)*X WE THEN SOLVE FOR THE EIGENVALUES OF THE X' = (M**0.5)*X. (SOME REFERENCES DESCRIBE THIS AS SOLVING FOR THE EIGENVALUES IN THE GENERALIZED COORDINATES. WE FINALLY SOLVE FOR X = (M**-0.5)*X'. IMPORTANT ASSUMPTIONS: K = IS SYMMETRIC AND REAL STIFFNESS MATRIX M = IS DIAGONAL MATRIX OF THE MASS USED IN CALCULATION. N-VECTOR OF REAL DOUBLE PRECISION VALUES. THE CALCULATED EIGENVALUES ARE PUT IN THIS VECTOR. N-VECTOR. WORKING VECTOR. X = N-BY-N MATRIX OF REAL DOUBLE PRECISION VALUES. THE CALCULATED EIGENVECTORS ARE PLACES COLUMN-BY-COLUMN IN THIS MATRIX. LAMBDA NDIM = THE ACTUAL DIMENSION OF THE ARRAYS IN THE MAIN PROGRAM. N = THE ORDER OF MATRICES K AND M. WANTX = LOGICAL. IF WANTX = .FALSE. THEN SYMEIG SKIPS THE CALCULATION OF THE EIGENVECTORS, SAVING SOME TIME. NDIM = 9 WANTX = .TRUE. N=9 READ MATRIX K, ROW-BY-ROW. ALSO, CALCULATE M*** *0. 5. IF YOU NEED TO WORK WITH THE MASS OF THE FLOOR PLUS PARTICIPANTS, YOU HAVE TO CHANGE THE MATRIX FMASS(I) IN THIS LOOP TO M(I) DO 250 I = 1, N READ(2, *) (K I, J), J=1, N) IM(I) = 1. 0D0 DSQRT(FMASS(I)) 250 CONTINUE CALCULATE (M**-0.5)*K DO 260 I = 1, N DO 260 J = 1, N PHI(I,J) = M(I)*K(I,J) SAVE THE STIFFNESS MATRIX IN THE ARRAY STIF(9,9) STIF(I,J)=K(I,J) 260 CONTINUE CALCULATE (M**-0.5)*K*(M**-0.5) DO 270 I = 1, N DO 270 J = 1, N K(I, J) = PHI(I, J)*M(J) 270 CONTINUE CALCULATE THE EIGENVALUES (AND EIGENVECTORS IF DESIRED). CALL SYMEIG (NDIM,N,K,LAMBDA,E,WANTX,PHI) SKIP THE EIGENVECTORS IF WE DIDN'T ASK FOR THEM IF (.NOT.WANTX) Go To 320 TRANSFORM THE EIGENVECTORS BACK FROM THE "GENERALIZED COORDINATES" DO 310 J = 1, N 0000 0000 000 0000 0000 0000000 310 320 330 340 360 350 370 390 380 410 400 395 154 D0 310 I = 1, N PHI(I,J) = PHI(I,J)*M(I) CONTINUE CONTINUE WRITE THE EIGENVECTORS INTO FILE='PHI'. THE EIGENVECTORS PHI(J), J=1,2,...,9, ARE PLACED ROWFBY-ROW IN THIS FILE. DO 330 J=1,9 ‘WRITE (6,340)(PHI(I,J),I=1,9) CONTINUE FORMAT(9E12.6) REWIND 6 DO 350 I=1,9 DO 360 J=1,9 MA(I,J)=TMASS (I,J) CONTINUE CONTINUE CALCULATE MNN(N)=[PHIT](N)*[MA]*[PHI](N),[PHIT](N) IS THE TRANSPOSITION OF NTH EIGENVECTOR. DO 370 I=1,9 READ(6,340)(PHIT(J),J=1,9) II=l JJ=9 KK=9 CALL MATRIX(PHIT,MA,MN,II,JJ,KK) II=1 JJ=9 KK=1 CALL MATRIX(MN,PHIT,MNN(I),II,JJ,KK) CONTINUE REWIND 6 CLEAR OUT THE MATRIX DA BEFORE STARTING DO 380 I=1,9 DO 390 J=1,9 DA(I,J)=0.0 CONTINUE CONTINUE SUM. OF (2*XIN*OMEGAN/MNNN)*[PHI]*[PHIT]) CALCULATE DA SUM. OF (2*XIN*OMEGAN/MNNN)*DAA) XIN=0.015 DO 395 I=1,9 READ(6,340)(PHIT(J),J=1,9) II=9 JJ=l KK=9 CALL MATRIX(PHIT,PHIT,DAA,II,JJ,KK) DO 4004.51,}1 9 D0 = DA(J,L)£DA(J,L)+(2*XIN*DSQRT(LAMBDA(I))/MNN(I))*DAA(J,L) CONTINUE CONTINUE CONTINUE REWIND 6 CALCULATE [DAMP] = [MA] * [DA] * [MA] II=9 JJ=9 KK=9 CALL MATRIX(MA,DA,DAM,II,JJ,KK) CALL MATRIX (DAM, MA, DAMP, II, JJ, KK) + ++ +++++4 +++++++++++++ 7. OTAL FORCES +++++++++++++ 155 CALCULATE [EA]=[M]* ANODE , FV = DAMP * VNOD ’ = * [FORCE]=[EAJ+[FV19[FX] 1 [ ] [ ] [ E] [FX] [K] [XNODE],AND DO 420 I=1,124 DO 430 J=1,9 XNODE(J)=DEF(I+2,J) VNODE(J)=VEL(I+2,J) ANODE(J)=ACCL(I+2,J) 430 CONTINUE II=9 JJ=9 KK=1 CALL MATRIX (MA, ANODE, FA, II, JJ, KK) CALL MATRIX(DAMP,VNODE,FV,II,JJ,KK) CALL MATRIX(STIF,XNODE,FX,II,JJ,KK) DO 440 J=1,9 FORCES(I,J)=FA(J)+FV(J)+FX(J) 440 CONTINUE WRITE(3,470) (FORCES(I,J),J=1,9) 420 CONTINUE DO 450 I=1,124 DO 460 J=1,9 TOTFORCE(I)=FORCES(I,J)+TOTFORCE(I) 460 CONTINUE WRITE(4,*)TOTFORCE(I) 450 CONTINUE 470 FORMAT(9F11.4) TO FIND THE MAX. FORCE ON THE FLOOR 000000 000 N=o FMAX=0.0 NMAX=0 DO 480 I=1,124 =N+1 IF (TOTEORCE(I).GT.EMAX) THEN FMAX=TOTFORCE(I) ( NMAX=N I GO TO 480 ‘ ENDIF ‘ 480 CONTINUE I WRITE(*,490)FNAME,NPART,DATE,STIME,ETIME,NODEBASE,NV,VMAX ( WRITE(S,490)FNAME,NPART,DATE,STIME,ETIME,NODEBASE,NV,VMAX WRITE(*,500)NODEBASE,MAX,FMAX,FMAX/NPART,FMAX/NPART/3.5 WRITE 5,500)NODEBASE,MAX,FMAX,FMAX/NPART,FMAX/NPART/3;5 490 FORMAT(l ,28X,'OUTPUT SUMM?RY},//,1X,'INPUT FILE NAME : ,A12,// *,1X,'NO. OF PART. : ',I , . *,1X,'DATE OF TEST : ',A12,//,1X,'TEST STARTED AT : ,A12,//,1X, *'TEST ENDED AT : ',A12/;,1X,'MAX. VOLT. (NODE ',I1,') IS SAMPLE *NO. ' I3 2X '=' 1X,E10.2, ) 500 EORMAT(1X,'MAX.'PEAK (NODE ',I1,') IS SAMPLE NO.',I3,3X,'= ',F9.2 *,2X,'lb',//,1X,'THE AVERAGE DYNAMIC HUMAN LOAD IS = ,F9.2,3X, *'1b/PERSON',//,1X,'THE DESIGN EQUIVELANT LOAD IS = ,1x, *F7.2,4X,'PSF',/) STOP END C C**************************************************************************** C SUBROUTINE SYMEIG (NDIM,N,A,D,E,WANTX,X) INTEGER NDIM,N LOGICAL WANTX REAL*8 A(NDIM,NDIM),D(NDIM),E(NDIM),X(NDIM,NDIM) REAL*8 ALPHA,BETA,GAMMA,KAPPA,AIJ,T,C,S,F,DABS,DSQRT S EIGENVALUES AND EIGENVECTORS OF REAL SYMMETRIC MATRICES FHOHUELINEAR MATHEMATICS FOR ENGINEERING APPLICATION" BY JOEL H. FERZIGER (1978). NDIM = DECLARED ROW DIMENSION OF A (AND X) N = ORDER OF.A = — - TRIX DOUBLE PRECISION) ‘A .A IS ASEUMED TO BE SYMMETRIC. ONLY THE DIAGONAL AND LOWER 0000000000 000000000000000000000000 000 0000 0000 WANTX = 156 TRIANGULAR PORTIONS OF A.ARE USED. IF POSSIBLE, PUT THE BIGGEST ELEMENTS IN THE UPPER LEFT CORNER. UPPER TRIANGLE PORTION OF A.IS UNALTERED UNLESS A.AND X ARE IN THE SAME N-DIMENSIONAL VECTOR, (DOUBLE PRECISION) STORES THE OUTPUT EIGENVALUES. N-DIMENSIONAL VECTOR, (DOUBLE PRECISION) USED FOR WORK SPACE. .A LOGICAL VARIABLE. IF WANTX = .TRUE. THEN THE EIGENVECTORS ARE CALCULATED. IF WANTX = .FALSE. THEN THE EIGENVECTORS ARE NOT CALCULATED. N-BY-N MATRIX (DOUBLE PRECISION) IF WANTX = .TRUE. THEN THE OUTPUT X(I,J) IS THE EIGENVECTOR ASSOCIATED WITH EIGENVALUE D(J). CALCULATES THE EIGENVECTORS AND REPLACES THE INPUT MATRIX, A, WITH ITS EIGENVECTORS. HOUSEHOLDER TRIDIAGONALIZATION D(1) = A(1,l) IF (N.LE.2) GO TO 8 NMl = N-l DO 7 K = 2, NMl FIND REFLECTION WHICH ZEROS A(I,K-1), I=K+1,...,N ALPHA = 0. DO 1 I = K N D(I) ='A(I,K—1) ALPHA = ALPHA + D(I)*D(I) CONTINUE ALPHA = DSQRT(ALPHA) IF (D(K).LT.0.) ALPHA = -ALPHA D(K) = D(K) + ALPHA BETA = ALPHA*D(K) IF(BETAJEQ.0.) GO TO 6 APPLY REFLECTION TO BOTH ROWS AND COLUMNS OF REST OF A USE AND COMPUTE ONLY LOWER TRIANGLE KAPPA = 0. DO 3 I = K, N GAMMA = 0. DO 2 J = K, N IF (I.GE.J) AIJ A(I,J) IE (I.LT.J) AIJ A(J,I) GAMMA = GAMMA + AIJ*D(J) CONTINUE E(I) = GAMMA/BETA KAPPA = KAPPA + D(I)*E(I) CONTINUE KAPPA = 0.5*KAPPA/BETA DO 5 I = K, N E(I = (I) - KAPPA*D(I) ) J - K, I DO 4A(I,J) = A(I,J) — D(I)*E(J) - E(I)*D(J) CONTINUE CONTINUE .A(K, D(K) K-1) = D(K) =.A(K,K) K) = BETA =.A(N,N) = A(N,N—1) ACCUMULATE TRANSFORMATIONS PRODUCES X SO THAT XT*(INPUT A)*X IS TRIDIAGONAL IF(. X(N, NOT.WANTX) GO TO 20 N) = 1. 11 12 13 15 0000 20 21 000 22 000 000 000 23 24 157 DO 15 KB = 1,NM1 K = N-KB BETA = A(K,K) DO 11 I = K, N X(I,K) = 0. X(K,I) = 0. CONTINUE X(K,K) = . IF (K.EQ.1 .OR.BETA, EQ.0.) GO TO 15 DO 14 J = K, N GAMMA = 0. DO 12 I = K, N GAMMA = GAMMA + X(I,J)*A(I,K-1) CONTINUE DO 13 I = K, N X(I,J) = X(I,J) - GAMMA*A(I,K-1) CONTINUE CONTINUE CONTINUE TRIDIAGONAL QR ALGORITHM IMPLICIT SHIFT FROM LOWER 2-BY-2 DO 27 MB = 2, N M = N+2-MB MM1 = M-l ITER = 0 L = 1 E(L) = 0. FIND L SUCH THAT E(L) IS NEGLIGIBLE L M S DABS(D(L-l)) + DABS(D(L)) T S + DABS(E(L)) IF (T.EQ.S) GO TO 23 L L-l IF (L.GE.2) GO TO 22 IF E(M) IS NEGLIGIBLE, THE D(M) IS AN EIGENVALUE, SO... IF (L.EQ.M) GO TO 27 IF (ITER.GE.30) GO TO 27 ITER = ITER + 1 FORM IMPLICIT SHIFT T (D(M—l) - D(M))/(E(M) + E(M)) S DSQRT(1. + T*T) IE (T.LT.0.) S = -S S = D(M) - E(M)/(T+S) E(L) = D(L) - S E = E(L+1) CHASE NONZERO F DOWN THE MATRIX DO 26 J = L,MM1 T = DABS(E(J)) + DABS}F) ALPHA = T*DSQRT((E(J) T)**2 + (F/T)**2) C = E J)/ALPHA = F ALPHA TA = S*(D(J+1) — D(J)) + 2.0*C*E(J+1) ) = ALPHA + ) = E(J+l) - C*BETA *BETA = D(J) + T ) = D(J+1) - T J.EQ.MM1) GO TO 24 E = S*E(J+2) E(J+2) = —C*E(J+2) S BE E(J E(J 1 T = S D(J) D(J+1 IF ( IF(.NOT.WANTX) Go TO 26 DO 25 I = 1, N T = X(I,J) X(I,J) = C*T + S*X(I,J+1) X(I,J+l) = S*T - C*X(I,J+l) 0000000 25 26 27 30 10 158 CONTINUE CONTINUE GO TO 21 CONTINUE RETURN END ****************************************************************** SUBROUTINE MATRIX.FOR ****************************************************************** SUBROUTINE MATRIX(A,B,C,M,N,L) INTEGER M, N, L DIMENSION A(M,N),B(N,L),C(M,L) DO 10 I=1,M DO 20 K=1,L C(I,K)=0.0 Do 30 J=1,N C(I,K)=C(I,K)+A(I,J)*B(J,K) CONTINUE CONTINUE CONTINUE RETURN END APPENDIX B LOAD SIMULATION PROGRAMS l .SIIHLSJWZR $STORAGE:2 O 0000000*********$**********X-li-li-X-fl-X-ifi-fl- 000 0000 $LARGE * *-**-**-**-** **-**-**-**-**-**-**-*x-**-* *********************** *********************** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** OCT 23 , 1992 BASSEM MICHIGAN STATE UNIVERSITY EAST LANSING, MI 48823 U.S.A ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** THIS PROGRAM READS THE PROCESSED *.CRL FILES AND PERFORMS THE FOLLOWING TASKS: K. KHAFAGI FOR.A CERTAIN ACTION *it-X-X-X-li-X-fl-X’N-i-Xvi-*IQ-i-fl-X-Ifl-Jffl-l-fl-i-i-i-i-JP ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. INPUT/OUTPUT MANAGEMENT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ PROGRAM hls.FOR INTEGER WEIGHT(1000) INTEGER WEIGHT(lO) CHARACTER INPUTS*40,FOUT*12 CHARACTER FNAME1*2,FNAME*12,0UTPUTS*40 COMMON /CMP/ M,N,FS,IWIN,L,NFFT,DF COMMON /FILES/ INPUTS,OUTPUTS,LINP,LOUT common statement to save the array u from function uni between calls common u(l7) in order to assure that no sequence of random numbers 9 enerated by function uni will be used twice, the seed array for uni (array u) 18 saved from run to run in file 'Seedarray'. 160 in order to get 161 c started, for the first run only, a separate program is run the only c purpose of which is to invoke the function rstart that is in uni. c rstart initializes the seed array u. for all subsequent runs the 0 seed array from the prev1ous run will be used instead of invoking c rstart by running the separate program. c open(l,file='sj-hls.dat') do 4 i=l,17 . read(l,9) u(1) 9 format(1x,e27.20) 4 continue C c DATA.INPUT PROVIDED BY THE USER c c 1000 format (/1X,a,' ? '$) PRINT *,' 1 I PRINT *,' 1 I C write (*,1000) 'NO. of Persons Performing this Actions' 5 read (*,*) NPART if (NPART.lt. 0.) then write (*.'(/a)') * ' *** ERROR *** NPART is less than 0 - Reenter' go to 5 else if (NPART.gt.1000.) then wrlte (*,'(/a) ) * ' **In this BETA.version...No. of Part. can not exceed 1000' go to 5 end if C c . . . . 16 write (*,1000) 'Do you Want to input Part1c1pants Weight..(Y/N) read (*,'(a)') IWT IF((IWT.EQ.'Y').or.(IWT.EQ.'y')) THEN C C I do 6 i3=l,npart . 239 write (*,1010) I] . . ' ' ' 1010 format (/1X,fWeight of Part1c1pant No.: I4, ? $) read (*,*) weight(ij) if (weight(13).le. 0.0) then write (*,'(la)') , , ' * ' *** ERROR *** Weight is less than 0 - Reenter goto 239. elseif(weight(11).gt.300.0) then wrlte (":(/a2') ' ' 300 lb R t ' * ' Sorry,!!! My limitation is a Weight less than . een er goto 39. elseif(weight(il).1t.30.0) then write ( ,'(/a)') . . "' | * ' Sorry,!!! Infants were not conSIdered 1n HLS ..... Reenter goto 239 6 endif continue else IF((IWT.EQ.'N').or.(IWT.EQ.'n')) THEN C C C DO 17 In=1,NPART c c Weight r.v.: : a normally distributed deviation from the nominal value Will c be used. C (mean=+l73 lbs, std.dev.=16 lbs) for men C (mean=+l43 lbs, std.dev.=20 lbs) for women 162 c see REFRENCE: ..... c rnum=rnor(0) AVGWT=160.0 STDWT=17.0 18 weight(in)=AVGWT+rnum*STDWT if((weight(in).lt.0.0).or.(weight(in).gt.300.0))go to 18 if(npart.le.15) print *,'Weight of Participant No.:',In,' is = ' l,weight(in),' lbs.’ 17 continue else goto 16 end if 000 fnamel='SJ' WRITE (FNAME,40) FNAMEl WRITE (FOUT,41) FNAMEl 40 FORMAT (A2,'-HLS.DAT') 41 FORMAT (A2,'-HLS.OUT') FLAG FOR THE OPEN SCREEN OFLAG=1 Read INPUT Data from PARAMETER FILE "hls.DAT" 000 000 READ (1,10) INPUTS READ (1,10) OUTPUTS READ (1,3) BW 3 format (f3.l) na=1 n=90 10 FORMAT (A40) DETERMINE THE LOCATION OF INPUT AND OUTPUT FILES 1.0 INPUTS 000000 LINP=LEN TRIM(INPUTS) IF (LINPTNE.0) THEN INPUTS=INPUTS(1:LINP)//'\' LINP=LINP+1 END IE 2.0 OUTPUTS 000 LOUT=LEN TRIM(OUTPUTS) IF (LOUTTNE.0) THEN OUTPUTS=OUTPUTS(1:LOUT)//'\' LOUT=LOUT+1 END IF OPEN NEEDED FILES OPEN (7,FILE=OUTPUTS(1:LOUT)//'hls.OUT') OPEN (4,FILE=OUTPUTS(1:LOUT)//fout) 00 call sjump(Npart,weight,fname,oflag,avgwt,stdwt) 000000 00 U 0 Z 163 l , lPRINT *,' SUCCESSFUL RUN !!!! HERE IS THE SAVED OUTPUT 1 I 1 I l I PRINT *,' TWO'FILES WERE GENERATED AS OUTPUT ...... l I 1 I PRINT *,' FILE NO. 1 ... NAME : SJ-HLS.OUT l lPRINT * ' ' CONTENT : FORCE-TIME HISTORIES l I 1 I PRINT *,' FILE NO. 2 ... NAME : HLS.OUT I PRINT *,' CONTENT : STATISTICAL OUTPUT 1 I 1 I I PRINT * ' HLS.FOR, V8.1.0, NOV. 1992 ...... B. KHAFAGI, MSU I END C END OF THE MAIN PROGRAM DATA—P1.FOR C####################################################################### C * * ~k * * * ~k * ~k * * * 'k * 'k * * 'k 'k * 'k * * * 'k * * 'k * * -k * * 'k 'k * * * * * * * * * * ~k * 'k * * * 'k * * * * 'k 'k 'k * * 'k * * * * * * * ~k * ~k * * ****iv):*********************************************************** * ++++++++++++++++++++++++++++++++++ ++++++++++++++++ SUBROUTINE SJ (SINGLE JUMP) +++++++++++++++ ++++++++++++++++++++++++++++++++++ 000+000 SUBROUTINE s'um (N art weight fname oflag,av wt,stdwt) DIMENSION avgtIB)? sgdt(8),avga(8),stda(8),avg%(5),stdf(5) DIMENSION allin(22,5),ALLOUT(22,5),d1ff(22,5) DIMENSION COVt(7I7)'COVIé358)ICEY§(3)5)Ih(?889)X(8) DIMENSION corf 5 crv cor cora , DIMENSION FXT(( 5 )FYT(90),'FZT(90), TIMET(90), AT(8), TT(8) DIMENSION A(10 T(1000,8),PHI(1000),ALL(22,1000),DATA(1000) 9 0 ) dimension Error 5;CONTROL(90),e(90) INTEGER WEIGHT( C DIMENSION A(lO, T(10,8),PHI(10),ALL(22,10),DATA(10) C INTEGER WEIGHT( LOGICAL*1 XABS 5 O 0 ( 1 8 1 I ,8 90 00 )r 0) 164 CHARACTER INPUTS*40,outputs*40 CHARACTER FNAME*12 COMMON /CMP/ M,N, FS, IWIN, L, NFFT, DF COMMON /FILES/ INPUTS, OUTPUTS, LINP, LOUT common statement to save the array u from function uni between calls common u(l7) Read Spectral estimation parameters for the error term var=0.0244720 frmin=0.0 frmax=10.0 FS=37. dtt=1.0/FS nfft=128 TIMESTEP= 1. O/FS do 23 i=1,n timet(i) =timet(i-1)+timestep if(i. eq. 1) timet(i)= continue 000 000 O) ++++++++++++++++++++++++++++++++++ 3. FILES' PROCESSING ++++++++++++++++++++++++++++++++++ ++++++++++++++++ ++++++++++++++++ 0000000N npt=7 npa=8 npf=5 . . . ' _ . read (1,110) avgph1, stdphi,(avgt(1),1=l,npt),(stdt(1),1=l,npt) do 112 i=1, npt read (1,114) (cort(i, j), j=—1,npt) continue . . read (1,116) (avga(i), i=1,npa),(stda(i),i=1,npa) do 118 i=—1,npa read (1,119) (cora(i, j), j= -1, npa) continue read (1,122) (avgf(i), i=—1,npf),(stdf(i), i=1, npf) do 124 i=1, npf read (1,126) (corf(i, j), j= -1, npf) t. SnichnirtWE(137.3;S/,7(f7.3),/,7(f7.3)) 3),/.8(f7.3)) format format 3)) 8333/,5(f7.3)) 8) ( ( format E ( ( 112 118 124 110 114 116 119 122 126 136 format format format HIU'IU'ICDQQ HAAAAA NHH‘hHH‘hI'h CDQQ\.'I\)Q whi hi V {MOO-ti: 167 168 169 allin(i+ allin(i+ continue do 168 i= allin(i+9, allin(i+9, continue do 169 i= 1) 2) 1) 2) :1, 2,2 l, 1 2 —1, allin(i+17,1)=avgf allin(i+17, 2)=stdf continue (1'4 'O'OS QAQ t t Ia-H- AA AA P-IJ- vv vv ) ) (i (i 165 C c c now the log normal varables are calculated c c ++++++++++++++++ hase lag ccc phivx2=(stdp i/avgph1)**2.0 ccc stdphi=sqrt(log(1.0+phivx2)) ccc avgphi=1og(avgphi)-( og(1.0-phivx2)/2.0) c ++++++++++++++++ loads as log normal ccc do 47 i=1,npa ccc if((i.e3.3).or.(i.eq.6)) then ccc vx2=(st a(i)/avga(1))**2.0 ccc avga(i) =log(avga(i))-(lo (1.0-vx2)/2.0) ccc stda(i) =sqrt(1og(1.0+vx2?) ccc endif 47 continue ccc do 45 i=1,npf . . ccc if((i.eq.1).or.(1.eq.3).or.(1.eq.5)) then ccc vx2=(stdf(i)/avgf(1))**2.0 ccc avgf(i) =log(avgf(1))-(log(1.O-vx2)/2.0) ccc stdf(i) =sqrt(log(l.0+vx2)) ccc endif 45 continue C I C now find the cov. matrix for amplit. (a), t1me t, force f . C by multiplying correlation coeff1c1ent times st.dev (1), (j) c c do 125 i=1,npt do 127 j=1,npt . . . ' covt(i,j)=cort(1,j)*stdt(1)*stdt(j) 127 continue 125 continue do 135 i=1,npa do 137 j=1,npa . . . . cova(i,j)=cora(1,j)*stda(1)*stda(j) 137 continue 135 continue do 155 i=1,npf do 157 j=1,npf . . . . covf(i,j)=corf(1,j)*stdf(1)*stdf(j) 157 continue 155 continue 0 c DO 380 K=1,Npart C ZERO THE FLAG Zi (USED IN THE LOAD MODELING PROCESS) C Zl=0.0 22:0.0 Z3=0.0 Z4=0.0 25:0.0 Z6=0.0 Z7=0.0 TIl=0.0 T12=0.0 now pick the random variables that represent the control points response lag phi a log-normally distributed deviation from the nominal value will b d. (ge:::+0.28 s, std.dev.=0.091bs) for men see REFRENCE: ..... 0000000000 IL CCC 0000 CCC CCC CCC CCC OOOOO 0000000000000 166 rnum=rnor(0) an=avgphi+rnum*stdphi avp=exp(an) phi(k)=an each of the time control will be assumed at this point noramllg distributed With correlated random variables called for time t -t7 call hmatrix (covt,crv,npt,npt,h,x) tt(l)=0.00 tt(2)=avgt(1)+crv(l) tt(3)=avgt(2)+crv(2) tt(4)=avgt(3)+crv(3) tt(5)=avgt(4)+crv(4) tt(6)=avgt(5)+crv(5) tt(7)=avgt(6)+crv(6) tt(8)=avgt(7)+crv(7) call hmatrix (cova,crv,npa,npa,h,x) at(1)=avga(1)+crv(l) at(2)=avga(2)+crv(2) at(3)=avga(3)+crv(3) at(3)=exp(av a(3)+crv(3)) at(4)=avga(4 +crv(4) at(5)=avga(5)+crv(5) at(6)=exp(avga(6)+crv(6)) at(6)=avga(6)+crv(6) at(7)=avga(7)+crv(7) at(8)=avga(8)+crv(8) call hmatrix (covf,crv,npf,npf,h,x) fxx=avgf(l)+crv(1) fxx=exp(avgf(1)+crv(l)) fxm=avgf(2)+crv(2) fyx=exp(avgf(3)+crv(3)) fyx=avgf(3)+crv(3) fym:avgf(4)+crv(4) DO 195 I=l,n CALCULATE THE VALUES OF THE CONTROL MODEL CONTROL(I) FOR AMPLITU ERROR BY USEIN THE SAMPLE TIME CONTROL POINTS INSTEAD OF MEAN TI CONTROL POINTS. IF (TIMET(I).LT.AVP) THEN CONTROL(I)=0.0 GO TO 195 ELSE IF (TIMET(I).LT.(TT(2)+AVP)) THEN STEP=((AT(2)-AT(1))/(TT(2)-TT(1)))/FS CONTROL(I)=CONTROL(I-1)+STEP 21:21+1. IF (Zl.EQ.l.0) CONTROL(I):AT(1) GO TO 195 ELSE IF (TIMET(I).LT.(TT(3)+AVP)) THEN AL USE A SECOND DEGREE PARABOLA TO MODEL THIS PART. THE INITI CONDITIONS FOR THE PARABOLA IS A.VALUE OF A(Z) AT TIME—0 (AT T2) AND A.VALUE OF.A(3) AT TIME=dt (AT T3) AND SLOPE OF 0 AT 0 F(t) =.ATA2 + BT + C AO=dA/dT“2 BO=0 CO=A2 dT=T3-§3 dA=A3- T(6) IS THE REFRENCE INITIAL TIME DT=TT(3)-TT(2) DA;AT(3)-AT(2) AO=DA/DT**2.0 BO=0.0 CO=AT(2) CC CC 000000000000 195 196 0 0000000 167 CONTROL(I)=AO*TI1**2+BO*TIl+CO TIl=TI1+1.0/FS STEP=((AT(3)-AT(2))/(TT(3)-TT(2)))/FS CONTROL(I)=CONTROL(I-1)+STEP zz=zz+1. IF (ZZ.EQ.1.0) CONTROL(I)=AT(2) GO TO 195 ELSE IF (TIMET(I).LT.(TT(4)+AVP)) THEN STEP=((AT(4)-AT(3))/(TT(4)-TT(3)))/FS CONTROL(I)=CONTROL(I—1)+STEP z3=z3+1. IF (Z3.EQ.1.0) CONTROL(I)=AT(3) GO TO 195 ELSE IF (TIMET(I).LT.(TT(5)+AVP)) THEN IF ((Z4.EQ.0.0).AND.(CONTROL(I-l).LT.-1.0)) CONTROL(I-1)=AT(4) STEP=((AT(5)—AT(4))/(TT(5)-TT(4)))/Fs CONTROL(I)=CONTROL(I-l)+STEP z4=z4+1. IF (Z4.EQ.1.0) CONTROL(I)=AT(4) GO TO 195 ELSE IF (TIMET(I).LT.(TT(6)+AVP)) THEN STEP=((AT(6)-AT(5))/(TT(6)-TT(5)))/FS CONTROL(I)=CONTROL(I—1)+STEP zs=zs+1. IF (ZS.EQ.1.0) CONTROL(I)=AT(5) GO To 195 ELSE IF (TIMET(I).LT.(TT(7)+AVP)) THEN USE A SECOND DEGREE PARABOLA TO MODEL THIS PART. THE INITIAL CONDITIONS FOR THE PARABOLA IS A.VALUE OF A(6) AT TIME=0 (AT T6) AND A VALUE OF.A(7) AT TIME=dt (AT T7) AND SLOPE OF 0 AT dT F(t) = ATAZ + BT + C AO=[A6-A7]/dT“2 BO=-2*A*dT CO=A6 dT=T7-T6 T(6) IS THE REFRENCE INITIAL TIME DT=TT(7)-TT(6) Ao=(AT(6)-AT(7))/DT**2.0 BO=-2*AO*DT co=AT(6) CONTROL(I)=AO*TI2**2+BO*TIZ+CO TIZ=T12+1.0/FS GO To 195 ELSE IF (TIMET(I).LT.(TT(8)+AVP)) THEN STEP=((AT(8)-AT(7))/(TT(8)-TT(7)))/FS CONTROL(I)=CONTROL(I—1)+STEP z7=z7+1. IF (Z7.EQ.1.0) CONTROL(I)=AT(7) GO To 195 END IF CONTROL(I)=0.0 CONTINUE do 196 i=l,8 t(k,i)=tt(1) a(k,i)=at(1) CONTINUE +++§++%#+%%%4+%:++4+%4++%+%%4+++++ ++++++++++++++++ 5. ERROR MODELENgLLLIILLL +++++++++++++++ +++++++++++++++++:L++++...11....T. call simu (VAR,FrMIN,FrMAX,DTt,NFFT,error) DO 210 I=1,n I; 0000000 000 0000000000 210 240 250 260 300 310 320 330 340 168 e(i)=error(i) FZt(I)=CONTROL(I)+e(i) Fxt(I)=e(i)/10.0 Fyt(I)=e(i)/30.0 ti=t1mestep*i-an cut =ti-tt(6) if((ti.ge.tt(4)).and.(ti.1e.(tt(5)+0.03)))then FZt(I)=CONTROL(I) endif 1f((cut.lt.timestep).and.(ti.ge.tt(6)))then fxt(i)=fxx fxt(i-1)=fxm EYE“??? yt 1- = endif ym CONTINUE +++++++++++++++++++++++++++++++++++++++ +++++++++++++ 5. START THE FINAL OUTPUT +++++++++++++ +++++++++++++++++++++++++++++++++++++++ PRINT TIME VS FORCES IN ACTUAL VERSUS ESTIMATED MODEL, AND ERROR IF (NPART.GT.10) THEN if(k.eq.1) then PRINT * ' FOR MORE THAN 10 PART., LOAD HISlTORIES ARE NOT SAVED' I PRINT *,' HOWEVER, CONTROL POINT STATISTICS ARE STILL SAVED.‘ PRINT *,' Press any key when ready.‘ read(*,*) endif ELSEIF (NPART.LE.20) THEN WRITE (4,240) . . FORMAT (///,10X,'TIME ',3X,'Z-FORCE Ratio ',2X,'MODEL Ratio',2X,' 1ERROR ARRAY'/) DO 250 I=1,90 WRITE (4,260) TIMEt(I),FZT(I),CONTROL(I),E(I) CONTINUE FORMAT (7X,F8.3,2X,F14.6,2X,F14.6,2X,F14.8) WRITE (4,300) WEIGHT(K) FORMAT (//,1x,'WEIGHT(S) :',I5,) FINALLY, WRITE A SUMMARY ABOUT THE FORCES. WRITE (4 310) fname FORMAT (50x,' FORCES INFORMATION ',///,13X,'FINAL REPORT FILE N 1AME = ',A12///,5X,'POINT',5X,'TIME',5X,'AMP.'/) DO 330 I=1,8 WRITE (4,320) I,T(K,I),A(K,I) FORMAT (1X,I3,3X,F8.3,3X,F14.6) CONTINUE WRITE (4,340) PHI(K) FORMAT (1X,'Response Lag (PHI) = ',F7.4) ENDIF PHIK=PHI(K) CALL PLOTXY SUBROUTINE TO PLOT THE TIME VERSUS THE FORCE IN THE X, Y, Z DIRECTIONS. NOW FIND THE MAXIMUM VALUE OF THE FORCE FX,FY,FZ, AND USE EACH FOR THE RANGE OF SCREEN PLOTS. 3 00000 0000 000 000 FZMIN=0.0 CALL CALL CALL CALL CALL CALL MAX MAX MAX MIN MIN MIN 169 (FXT,FXMAX,TX,TIMESTEP) (FYT,FYMAX,TY,TIMESTEP) (FZT,FZMAX,TZ,TIMESTEP) (FXT,FXMIN,TXM,TIMESTEP) (FYT,FYMIN,TYM,TIMESTEP) (FZT,FZMIN,TZM,TIMESTEP) FWEIGHT=WEIGHT(K) IF (NPART. GT. 10) THEN IF (K. EQ. 1) PRINT * ,'NO OUTPUTS TO THE SCREEN .... 1LIMITATIONS . "' ' PRINT *, ELSE SIMULATION No. : ',K IF (OFLAG.GT.0.0) THEN CALL OPEN (OFLAG,BW) ELSE GO TO 350 END IF 50 CALL PLOTXY (TIMET, FZT, CONTROL, FZMAX, FWEIGHT, BW, K, Npart, FZMIN 1, 90, FNAME, AVP, TT, AT, 8, FxT, foIN, foAX, fyt, fymin,fymax) ENDIF +++++++++++ 8. Reporting control points to fIiE.T§i?LQETTILTTTI+ IF (K. EQ.1) WRITE (7,360) WRITE (7,370) K, WEIGHT(K), PHIK,(TT(J), J=2, 8),(AT(I),I= 11, 8),FXMAX, FXMIN, FYMAX, FYMIN, FZMAX 360 FORMAT (//, 25X,'OUTPUT SUMMARY', //, 73X,'No WT PHI T2 1 T3 T4 T5 T6 T8 A1 . A2 13 A4 A5 A6 A7T A8 FXmax FXmin FYmax l FYmin FZmax' ) 370 FORMAT (1X, I3, 1X,i4, 21(2X, F6. 3)) PREPARE DATA FOR THE STATISTICS PACKAGE .ALL(1, K)= float(WEIGHT(K)) ALL(6E K) ;:tT ALL(16,K)= ALL(17,K)=a ) ) ) ) ) ) ) 1) 2) 3) 4) 5) 6 7 8 ) ) ) ALL(18, K) =FXMax ALL(19, K) =FXMin ALL(20, K) =FYMax ALL(21, K) =FYMin ALL(22, K)=FZMAX DONE WITH ALL FILE PROCESSING 380 CONTINUE add counter for the no. of time to call mystat nms=22 XABS=.false. C C CALCULATE AND REPORT IMPORTANT STATISTICS FOR THE OBSERVED ACTION do 600 i=1,nms do 610 k=l,NPART data(k)=all(i,k) 610 continue if(npart.eq.l)go to 611 CALL MYSTAT (data,NPART,DMIN,DMAX,DMEAN,DSTD,DSE,XABS) 611 continue c c same estimated variables in a new arrat ALLOUT in which the parameter . . . . c l is the variable being processed (eg. weight, ph1,..) and the parameter . . c nout is the output variable number defined as follows: c nest=2 allout(i,1) =Dmean allout(1,2) =Dstd 600 continue C C REPORT THESE VALUES TO THE SUMMERY OUTPUT FILE C c c BLANK LINE T0 SEPARATE DATA c do 666 i=1,22 . . diff(i,1)=abs((allin(i,1)—ALLOUT(I 1))/a111n(i,l)*100.0) diff(i,2)=abs((allin(i,2)-ALLOUT(I,2))/allln(1,2)*100.0) 666 continue WRITE (7,633) 633 FORMAT (1X/) WRITE (7,361) 361 FORMAT 25X 'OUTPUT SUMMARY',//,3X,'No WT PHI T2 1 T3 (/T& ' T5 T6 T7 T8 A1 A2 13 A4 .A5 A6 .A7 .A8') WRITE (7,653) (ALLin(I,1),I=1,17) WRITE (7,650) (ALLOUT(I,1),I=1,17) WRITE (7,655) (d1ff(1,1),I=l,l7) WRITE (7,654) (ALLln(I,2),I=1,17) WRITE (7,650) (ALLOUT(I,2),I=1,17) WRITE (7,655) (d1ff(1,2),I=1,17) 650 FORMAT (' Sim',f6.0,f8.3,15( 8.3)) 653 FORMAT (1x,'Mean valuesf,/,' Act',f6.0,f8.3,15(F8.3)) 654 FORMAT (1x,'Std. Deviat10ns',/,' Act',f6.0,f8.3,15(F§.3)) 655 FORMAT (' Dif',3x,f3.0,'% ',f5.l,'% ',15(F5.l, % )/) C CALL ENDGRAPHICS () RETURN *ENE‘kiri'**************************** ‘- X- X- * )I- *- 1- 3(- * *- #- ***‘k‘k *1-*#—*s—»»~*x~*x-+ x- x- a(- * i- X- x. » g. )1- fl- * * 3‘. * * a(- *- * 2(- (70 »»-**-*x-*x-*x-*x-* $i-*fi-*fl-*i-*i-*fi-* IL 00000 00 ***-***lt-X-X-i-i-OO * 9: 'k *- * * * * * * * l. 4. + ++++++++++++++++ 1 . DATA ++++++++++++++++ % + SUBROUTINE OPEN (OFLAG,BW) LOGICAL*1 INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lPLOTMETHOD,INVERTX,INVERTY,TITLEONLY,SCRONLY,Y2PLOT,Y3PLOT INTEGER*2 BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, lTITLECOLOR,NUMBER,PLOTZCOLOR,FONTHEIGHT,FONTWIDTH,PLOT3COLOR CHARACTER TITLE*280,XAXISTITLE*80,YAXISTITLE*80 INCLUDE 'PLOT.DIF' INITIALISE THINGS FOR THE FIRST PLOT INITGRAPHICS=.TRUE. INITBORDER=.TRUE. INITWINDOW=.TRUE. INITAXES=.FALSE. INITFONTS=.TRUE. PLOTMETHOD=.FALSE. Y2PLOT=.FALSE. Y3PLOT=.FALSE. SCRONLY=.TRUE. FONTHEIGHT=28 FONTWIDTH=16 SET UP THE FANCY COLORS BORDERCOLOR=5 WINDOWCOLOR=7 LABELCOLOR=4 SET UP THE PLOT LIMITS XMIN=0.0 XMAX=1.0 YMIN=0.0 YMAX=1.0 POSITION THIS PLOT AT THE TOP LEFT OF THE SCREEN XBOTTOMLEFT=0.0 YBOTTOMLEFT=0.0 XTOPRIGHT=1.8 $¥$E§i§HT 1 ** s J - H L s . F o R **' IF(BW.EQ.1) THEN BORDERCOLOR=O WINDOWCOLOR=15 LABELCOLOR=1 END IF CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, Y2DATA NUMBER PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX, $§MINé§M£X?XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT, 3YTOPRIGHT,PLOTZCOLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR, 4PLOTCOLOR,TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY, 5FONTHEIGHT,FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,X3,Y3,NUM3) WAIT FOR USER TO HIT RETURN READ(*r") CALL ENDGRAPHICS() OFLAG=O RETURN END END OF MODULE SUBROUTINE OPENING ############ ###§#§#f#§#f#f#§#g#g#f#f#f#f#f#f#f#f#f#f#f#f#f#f#f#f#f#f#§#. . . . . . ********************************* SUBROUTINE PLOTXY ****************~k ************* ** **************** **************** »» *1—*x-*x—**-*» APRIL 15,1992 PLOTTING FORCE-TIME HISTORIES * ******i-************************* 172 *******'k*************************** * ****************in):*****-kink**************************************** * 000000 ’1' * 000 00000 000 ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA.INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SUBROUTINE PLOTXY (TIME,FZ,PR,FZMAX,WEIGHT,EW,K,NF,FZMIN,N, 1FNAME,PHIK,X3,Y3,NUM3,fx,foIN,foAX,fy,fymin,fymax) DIMENSION TIME(90), FZ(90), PR(90) DIMENSION X3(NUM3), Y3(NUM3),XNULL(8),YNULL(8) INCLUDE 'PLOT.DIF' LOGICAL*1 INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lPLOTMETHOD,INVERTX,INVERTY,TITLEONLY,SCRONLY,Y2PLOT,Y3PLOT INTEGER*2 BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, lTITLECOLOR,NUMBER,PLOT2COLOR,FONTHEIGHT,FONTWIDTH,PLOT3COLOR CHARACTER TITLE*280,XAXISTITLE*80,YAXISTITLE*80, lFNAME*12,TITLE1*80,TITLE2*80,TITLE3*80,TITLE4* 280,TITLE5*80,TITLE6*80,TITLE7*80,TITLE8*80,TITLE9*80,TITLE10*80, 3TITLE11*80,TITLE12*80,TITLE13*80 NUMBER=N INITIALISE THINGS FOR THE FIRST PLOT INITGRAPHICS=.TRUE. INITBORDER=.TRUE. INITWINDOW=.TRUE. INITAXES=.TRUE. INITFONTS=.TRUE. PLOTMETHOD=.TRUE. TITLEONLY=.FALSE. SCRONLY=.FALSE. Y2PLOT=.false. Y3PLOT=.TRUE. FONTHEIGHT=10 FONTWIDTH=5 SET UP THE FANCY COLORS DEPENDING ON THE VALUE OF CFLAG (0 OR 1) BORDERCOLOR=14 IF (CFLAG.EQ.0.0) THEN WINDOWCOLOR=8 CFLAG=1.0 ELSE WINDOWCOLOR=1 CFLAG=0.0 END IF AXISCOLOR=10 LABELCOLOR=11 PLOTCOLOR=14 PLOT2COLOR=4 PLOT3COLOR=15 TITLECOLOR=15 ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. FORCE IN THE x-DIR ++++++++++++++++ ++++++++++++++++++++++++++++++++++ XMIN=0.0 XMAX=2.5 XTICSPACE=O.5 CREATE THE ZERO AXIS USING THE ARRAY XNULL XNULL(8)=xmax YNULL(8)=0.0 _ . . Call scale subroutine to set Ymin, Ymax, Yticspace,XMIN,XMAX,... IL 000 000 000000 000 173 CALL SCALE (1.3*fxmin,l.3*fxmax,YMIN,YMAX,YTICSPACE) POSITION THIS PLOT XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.45 XTOPRIGHT= 0. 50 YTOPRIGHT= O. 90 GI ITA TITLE TITLE= Fx AS A RATIO TO WEIGHT ' XAXISTITLE=' T I M E ( s e C . )' YAXISTITLE='FORCE RATIO' USE THE BLACK AND WHITE COMBINATION IF SCREEN CAPTURING IS NEEDED IF(BW.EQ.1) THEN BORDERCOLOR=0 WINDOWCOLOR=15 AXISCOLOR=1 PLOT3COLOR=2 TITLECOLOR=1 ND IF READY TO CALL PROGRAM PLOT.FOR CALL PLOT (INITGRAPHICS, INITBORDER, INITWINDOW, INITAXES, INITFONTS, lTIME, Fx, PR, NUMBER, PLOTMETHOD, INVERTX, INVERTY, XMIN, XMAX, YMIN ,YMAX, 2XTICSPACE, YTICSPACE, XBOTTOMLEFT, YBOTTOMLEFT, XTOPRIGHT, YTOPRIGHT, 3PLOT2COLOR, BORDERCOLOR, WINDOWCOLOR, AXISCOLOR,LABELCOLOR, PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,XNULL,YNULL,NUM3) ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. RCE IN THE y—DIR ++++++++++++++++ ++++++++++++++++++++++++++++++++++ CALL SCALE (l.2*fymin,1.2*fymaX,YMIN,YMAX,YTICSPACE) POSITION THIS PLOT INITGRAPHICS=.FALSE. XBOTTOMLEFT=O.5O YBOTTOMLEFT=0.45 XTOPRIGHT=1.00 YTOPRIGHT=O.9O GIVE IT A TITLE TITLE=' Fy AS A RATIO TO WEIGHT ' XAXIS MT E= ' TIME ( s e c . YAXISTITLE= 'FORCEI RATIO' USE THE BLACK AND WHITE COMBINATION IF SCREEN CAPTURING IS NEEDED IF(BW.EQ.1) THEN BORDERCOLOR=0 WINDOWCOLOR=15 AXISCOLOR=1 LABELCOLOR=4 PLOTCOLOR=0 PLOT2COLOR=1 PLOT3COLOR=2 TITLECOLOR=1 END IF READY TO CALL PROGRAM PLOT.FOR I; 000 0000000 000 000 00000 000 40 50 174 CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, 1TIME,Fy,PR,NUMBER,PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAX, 2XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,XNULL,YNULL,NUM3) ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. FORCE IN THE z-DIR ++++++++++++++++ ++++++++++++++++++++++++++++++++++ Call scale subroutine to set Ymin, Ymax, Yticspace,XMIN,XMAX,... CALL SCALE (l.2*fzmin,l.2*fzmax,YMIN,YMAX,YTICSPACE) POSITION THIS PLOT XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.00 XTOPRIGHT=0.50 YTOPRIGHT=O.45 GIVE IT A TITLE TITLE=' FZ AS‘A RATIO TO WEIGHT ' XAXISTITLE=' T I M E ( s e C . )' YAXISTITLE='FORCE RATIO' USE THE BLACK AND WHITE COMBINATION IF SCREEN CAPTURING IS NEEDED IF(BW.EQ.1) THEN BORDERCOLOR=O WINDOWCOLOR=15 AXISCOLOR=1 LABELCOLOR=4 PLOTCOLOR=0 PLOT2COLOR=1 PLOT3COLOR=2 TITLECOLOR=1 END IF READY TO CALL PROGRAM PLOT.FOR CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lTIME,FZ,PR,NUMBER,PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAX, ZXTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR1AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,XNULL,YNULL,NUM3) ++++++++++++++++++++++++++++++++++ +++++++++++++++ 2. OUTPUT VALUES ++++++++++++++++ ++++++++++++++++++++++++++++++++++ just plot the title on this one XBOTTOMLEFT=O.50 YBOTTOMLEFT=0.00 XTOPRIGHT=1.00 YTOPRIGHT=O.45 INITGRAPHICS=.EALSE. TITLE=' CONTROL POINTS ' WRITE (TITLE1,40) FORMAT ('Single Jump With One Person') WRITE (TITLE2,50) WEIGHT FORMAT ('Weight= ',F4.0,' lbs') Find the max. value of force in the x,y directions (absolute force WRITE (TITLE3,60) K,NF I=0 WRITE (TITLE13,80) PHIK I; 175 1:1 WRITE (TITLE4,70) I,X3(1),Y3(1) I=I+ WRITE (TITLE5,70) I,X3(2),Y3(2) I=I+ WRITE (TITLE6,70) I,X3(3),Y3(3) =I+ WRIT? (TITLE7,70) I,X3(4),Y3(4) I= + WRITE (TITLE8,70) I,X3(5),Y3(5) = + WRIT? (TITLE9,70) I,X3(6),Y3(6) = + WRITE (TITLE10,70) I,X3(7),Y3(7) = + WRITE (TITLE11,70) I,X3(8),Y3(8) 6O FORMAT (IZ,' OF'I3,1X,'POINT TIME AMP.') 70 FORMAT (8X,IZ,2X,F8.4,2X,F10.5) 80 FORMAT (1X,'Response Lag (PHI) = ',F7.4) WRITE (TITLE12,90) 9O FORMAT ('HIT ENTER TO CONTINUE....') CALL OUTPUTS (XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT,TITLE, lTITLEl,TITLE2,TITLE3,TITLE4,TITLE5,TITLE6,TITLE7,TITLE8,TITLE9, 2TITLE10,TITLE11,TITLE12,TITLE13,CFLAG,BW) C C C ++++++++++++++++++++++++++++++++++ C ++++++++++++++++ 3. Top Title on the Screen ++++++++++++++++ C ++++++++++++++++++++++++++++++++++ C INITGRAPHICS=.FALSE. TITLEONLY=.TRUE. INITAXES=.FALSE. INITFONTS=.TRUE. FONTHEIGHT=28 FONTWIDTH=16 WINDOWCOLOR=7 TITLECOLOR=4 XBOTTOMLEFT=0.00 YBOTTOMLEFT=0.9O XTOPRIGHT=1.00 YTOPRIGHT=1.00 C TITLE='S I N G L E J U M P' E USE THE BLACK AND WHITE COMBINATION IF SCREEN CAPTURING IS NEEDED IF(BW.EQ.1) THEN WINDOWCOLOR=15 TITLECOLOR=1 END IF CALL PLOT (INITGRAPHICS,INITBORDER,INITWINDOW,INITAXES,INITFONTS, lTIME,FZ,PR,NUMBER,PLOTMETHOD,INVERTX,INVERTY,XMIN,XMAX,YMIN,YMAX, 2XTICSPACE,YTICSPACE,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT,YTOPRIGHT, 3PLOT2COLOR,BORDERCOLOR,WINDOWCOLOR,AXISCOLOR,LABELCOLOR,PLOTCOLOR, 4TITLECOLOR,XAXISTITLE,YAXISTITLE,TITLE,TITLEONLY,FONTHEIGHT, C 5FONTWIDTH,SCRONLY,Y2PLOT,Y3PLOT,PLOT3COLOR,XNULL,YNULL,NUM3) C WAIT FOR USER TO HIT RETURN ( only for testin the program g , otherwise let the program run Wlth no stops? READ (*,*) IF (K.EQ.NF+1) CALL ENDGRAPHICS () RETURN END END OF SUBROUTINE PLOTXY ################################################################### * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *****3"**************************** * **»»n() m: =23: *xq: 000 000000It-i-1I-1t-x-It-1-1-1-fl- 176 SUBROUTINE MAX MAY 05,1992 FINDING THE MAX. VALUE OF DATA * 'k * 'k 'k * * * 'k 'k * * * 'k 'k * 'k * 'k 'k ‘k * * * * 'k * * * * * * ~k * * 'k * * 'k 'k * 'k 'k 'k 'k * 'k * * * 'k * 'k 'k 'k 'k 'k 'k * * 'k -k * * 'k 'k *****-****‘k******************************************************** *x-**-**-**-** *1—*fl-*fi-+i-** ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SUBROUTINE MAX (DATA,FMAX,T,TIMESTEP) TO FIND THE MAX. FORCE DIMENSION DATAX90) T=0.0 N=0 FMAX=-1El4 NMAX=0 DO 10 I=1,9O N=N+1 IF (DATA(I).GT.FMAX) THEN FMAX=DATA(I) NMAX=N GO TO 10 END IF lO CONTINUE DOC) nooonnx-x-x-x-x-x-x-x-x-x-x-x-zt-OO # ################################################################## ********* ******************** ******************** T=REAL(NMAX)*TIMESTEP RETURN END END OF SUBROUTINE MAX ** *3? ************* SUBROUTINE MIN 1* * * * * ****it***************************** * * 9: *- * * * **-**-»x-+x-*:-»x-*§: =11: ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ SUBROUTINE MIN (DATA,FMIN,T,TIMESTEP) TO FIND THE MIN. FORCE DIMENSION DATA(90) T=0.0 N=0 FMIN=0.0 NMAX=0 DO 10 I=1,90 N=N+1 IF (DATA(I).LT.FMIN) THEN I; 10 o********+**»000 0 00000 ***$****#***#3¥2 10 20 177 FMIN=DATA(I) NMAX=N GO TO 10 END IF CONTINUE T=REAL(NMAX)*TIMESTEP RETURN END END OF SUBROUTINE MIN ################################################################### ################################################################### ********************************* ********************************* * SUBROUTINE outputs * * **************************~k****** *********9:*********************** * MAY 15,1992 Plotting the final values of Forces * * ********************************* ********************************* START OF MODULE outputs.for INCLUDE 'FGRAPH. FI' SUBROUTINE OUTPUTS (XBOTTOMLEFT, YBOTTOMLEFT, XTOPRIGHT, YTOPRIGHT, 1TITLE,TITLE1, TITLEZ, TITLE3,TITLE4,TITLE5,TITLE6,TITLE7,TITLE8, ZTITLE9, TITLEIO, TITLE11,TITLE12,TITLE13, CFLAG, BW) INCLUDE 'FGRAPH. FD' INTEGER*2 DUMMY,XWIDTH, YHEIGHT, 1X1, IYl, IX2, 1Y2, IX, IY, FONTHEIGHT, 1FONTWIDTH REAL*4 XW1,YW1,XW2, YW2,XBOTTOMLEFT,YBOTTOMLEFT,XTOPRIGHT, lYTOPRIGHT, YT CHARACTER TITLE*80, FONTCOMMAND*10, FONTSTRING*7, FONTFILE*12, TITLE1* 180, TITLE2*80, TITLE3*80, TITLE4*80, TITLE5*80, TITLE6*80, TITLE7*80, 2TITLEB*80, TITLE9*80, TITLE10*80, TITLE11*80, TITLE12*80, TITLE13*80 RECORD / VIDEOCONFIG / SCREEN RECORD / WXYCOORD / WXY RECORD / XYCOORD / POSITION COMMON SCREEN SET UP FONT FOR USE FONTFILE='HELVB.FON' FONTCOMMAND="T'HELV'" IF (REGISTERFONTS(FONTFILE).LT.0) THEN WRITE (*,10) FONTFILE FORMAT (' FONT FILE ',A12,' NOT FOUND IN PLOT') STOP END IF FONTHEIGHT=28 FONTWIDTH=16 WRITE (FONTSTRING, 20) FONTHEIGHT, FONTWIDTH FORMAT (' H', IZ. 2,'W', I2. 2, 'B' DUMMY=SETFONT(FONTCOMMAND//FONTSTRING) GET SCREEN DIMENSIONS XWIDTH=SCREEN.NUMXPIXELS YHEIGHT=SCREEN.NUMYPIXELS SET A.VIEW PORT, WITH CORNER COORDINATES DEFINED AS.A FRACTION OF THE TOTAL SCREEN SIZE (0.0,0.0) IS BOTTOM LEFT OF SCREEN (1.0,1.0) IS TOP RIGHT OF SCREEN TRANSLATE CORNER COORDINATES TO PIXEL COORDINATES IX1=INT(XBOTTOMLEFT*FLOAT(XWIDTH)) IX2=INT(XTOPRIGHT*FLOAT(XWIDTH)) IY1=YHEIGHT-INT(YTOPRIGHT*FLOAT(YHEIGHT)) IY2=YHEIGHT-INT(YBOTTOMLEFT*FLOAT(YHEIGHT)) CALL SETVIEWPORT (IX1,IY1,IX2,IY2) CREATE.A REAL COORDINATE WINDOW XW1=0.0 YW1=0 00000 000 000 0000 00000 178 XW2=1.0 YW2=1.0 DUMMY=SETWINDOW(.TRUE.,XW1,YW1,XW2,YW2) SET THE COLOR OF THE BACKGROUND WINDOW DEPENDING ON THE FLAG "CFLAG" ... NOTE THE SETUP HAS TO BE OPPOSITE TO THE PLOTXY IF (CFLAG.EQ.1.0) THEN DUMMY=SETCOLOR(8) ELSE DUMMY=SETCOLOR(1) END IF IF(BW.EQ.1) DUMMY =SETCOLOR(15) DUMMY=RECTANGLE W($GFILLINTERIOR,XW1,YW1,XW2,YW2) DRAW A BORDERIAROUND THE VIEWPORT DUMMY=SETCOLOR(14) IF(BW.EQ.1) DUMMY =SETCOLOR(0) DUMMY=RECTANGLE_W($GBORDER,XW1,YW1,XW2,YW2) This is the screen design of the data output TITLELENGTH=50 POSITION THE TITLE CENTERED (MORE OR LESS) ABOVE GRAPH YT=0.9*(YW2+YW1) CALL MOVETO W (0.9* (XW2+XW1) ,YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-(TITLELENGTH*FONTWIDTH)/2 IY=IY-FONTHEIGHT/2 enter the program.name (Title: from the main program) DUMMY=SETCOLOR(15) IF(BW.EQ.1) DUMMY =SETCOLOR(O) CALL MOVETO (IX,IY,POSITION) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) DUMMY=SETCOLOR(4) IF(BW.EQ.1) DUMMY =SETCOLOR(0) CALL MOVETO (IX,IY+2,POSITION) CALL OUTGTEXT (TITLE(1:TITLELENGTH)) CHANGE THE FONT TO A FIXED FONT TO CREATE THE FORCES' TABLE HERE, COURIER FONTS WILL BE USED. FONTFILE='courb.fon' FONTCOMMAND="T’courier'" FONTHEIGHT=12 FONTWIDTH=9 IF (REGISTERFONTS(FONTFILE).LT.O) THEN WRITE (*,10) FONTFILE STOP END IF WRITE (FONTSTRING,20) FONTHEIGHT,FONTWIDTH DUMMY=SETFONT(FONTCOMMAND//FONTSTRING) YT=0.8*(YW2+YW1) CALL MOVETO W (0.73*(XW2+XW1),YT,WXY) CALL GETCURRENTPOSITION (POSITION) IX=POSITION.XCOORD IY=POSITION.YCOORD IX=IX-(TITLELENGTH*FONTWIDTH)/2 IY=IY-FONTHEIGHT/2 I; 179 DUMMY=SETCOLOR(O) IF(BW.EQ.1) DUMMY =SETCOLOR(1) CALL MOVETO (IX,IY+10,POSITION) CALL OUTGTEXT (TITLE1(1:TITLELENGTH)) CALL MOVETO (IX,IY+25,POSITION) CALL OUTGTEXT (TITLE2(1:TITLELENGTH)) DUMMY=SETCOLOR(11) IF(BW.EQ.1) GOTO 22 CALL MOVETO (IX,IY+ll,POSITION) CALL OUTGTEXT (TITLE1(1:TITLELENGTH)) CALL MOVETO (IX,IY+26,POSITION) CALL OUTGTEXT (TITLE2(1:TITLELENGTH)) 22 CONTINUE FONTHEIGHT=10 FONTWIDTH=8 WRITE (FONTSTRING,20) FONTHEIGHT,FONTWIDTH DUMMY=SETFONT(FONTCOMMAND//FONTSTRING) DUMMY=SETCOLOR(14) IF(BW.EQ.1) DUMMY =SETCOLOR(0) CALL MOVETO (IX,IY+40,POSITION) CALL OUTGTEXT (TITLE3(1:TITLELENGTH)) DUMMY=SETCOLOR(15) IF(BW.EQ.1) DUMMY =SETCOLOR(1) CALL MOVETO (IX,IY+55,POSITION) CALL OUTGTEXT (TITLE4(1:TITLELENGTH)) CALL MOVETO (IX,IY+68,POSITION) CALL OUTGTEXT (TITLE5(1:TITLELENGTH)) CALL MOVETO (IX,IY+81,POSITION) CALL OUTGTEXT (TITLE6(1:TITLELENGTH)) CALL MOVETO (IX,IY+94,POSITION) CALL OUTGTEXT (TITLE7(1:TITLELENGTH)) CALL MOVETO (IX,IY+107,POSITION) CALL OUTGTEXT (TITLE8(1:TITLELENGTH)) CALL MOVETO (IX,IY+120,POSITION) CALL OUTGTEXT (TITLE9(1:TITLELENGTH)) CALL MOVETO (IX,IY+133,POSITION) CALL OUTGTEXT (TITLE10(1:TITLELENGTH)) CALL MOVETO (IX,IY+146,POSITION) CALL OUTGTEXT (TITLE11(1:TITLELENGTH)) DUMMY=SETCOLOR(2) IF(BW.EQ.1) DUMMY =SETCOLOR(O) CALL MOVETO (IX,IY+164,POSITION) CALL OUTGTEXT (TITLE13(1:TITLELENGTH)) DUMMY=SETCOLOR(O) CALL MOVETO (IX,IY+185,POSITION) CALL OUTGTEXT (TITLE12(1:TITLELENGTH)) DUMMY=SETCOLOR(4) CALL MOVETO (IX,IY+186,POSITION) CALL OUTGTEXT (TITLE12(1:TITLELENGTH)) ALL DONE CALL UNREGISTERFONTS () RETURN END 000 ****************** ****************** w+ ** u-x- *» ** * * *$ *$ SUBROUTINE SCALE ********************** ***~k**************~k*** n~ a» x- 1. 2+2:- *3!- 2+:- x-x- »* $a+ *$ ***-***fi-fl-X- *****$*** :9 x- x- a(- AUG 22,1992 VERTICAL PLOT SCALE I; 180 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ***************************************************************** * ++++++++++++++++++++++++++++++++++ ++++++++++++++++ 1. DATA INPUT/OUTPUT ++++++++++++++++ ++++++++++++++++++++++++++++++++++ *-*fl-* ()0(3C)* *X-* SUBROUTINE SCALE (FMIN,FMAX,ST,FIN,TICK) This routine determines the starting and ending values to be used in a full scale plot. . . OUTPUT : ST = Starting value; FIN = Ending value; TICK . THIS IS A MODIFIED VERSION OF.A PROGRAM WRITTEN ORIGINALLY BY R. HARIC REAL*8 U,C TT=ALOG10(FMAX-FMIN) MT=IRINT(TT) TT=10.**(TT-MT) IF (TT.LE.2) THEN TT=2. ELSE IF (TT.GT.5) THEN TT=10. ELSE TT=5. END IF TICK=TT*10.**(MT-l) U=l.00001*TICK C=DLOG10(U) IF (C.LT.0) THEN L=IDINT(C)-1 ELSE L=IDINT(C) END IF C=DBLE(10.)**L TICK=C*IDNINT(U/C) N1=IRINT(FMIN/TICK) N2=IRINT(FMAX/TICK)+1 ST=N1*TICK FIN=N2*TICK RETURN END INTEGER FUNCTION IRINT(X) IF (X.LT.O) THEN IRINT=INT(X)-l ELSE IRINT=INT(X) END IF RETURN END ############################### ############################### function uni(jdum): uses mcgill ufierduper to generate uniform random numbers between 0 and 1. t e argument (jdum) is just a dummy argument - i.e. it isn't used in the subpro ram. for more information contact dr. george marsaglia was ington state university pullman,washington. 0000000 #flm mum 2 ##########################§# # # ######################### figgfifigfi s the function rstart is invoked once at the beginning of the main program to initialize uni. the arguments for rstart can be chosen arbitrarily. ******************************************************************** function uni(jdum) . integer nbits,i,j,m,11 common u(l7) lg; 000 0000000000 00 181 data nbits,i,j,m,i1/48,17,5,32707,1971/ . THIS SEED WAS USED TO SIMULATE random motions jdum:0 uni=u(i)-u(j) _ if(uni.lt.0.) uni=uni+l. u(i)=uni i=i-1 ifgi.eq.0) i=1? if(j.eq.0) j=17 return rstart initializes uni entry rstart(is,js) ' _ . . initializes u(1) to u(l7) bit-by-bit, uSing fibonaccl mod m, lag 2 1 =13 i3='s do 2 I=1,17 s=.5 u(l)=0. do 2 ll=1,nbits k=i3-il il=i2 i2=i3 i3=k if(i3.gt.0) go to 2 i3=i3+m u(l)=u(l)+s s=.5*s rstart=is*10000+js return end ********************************************************************* function rnor(j): generates normally distributed random numbers be- tween -infinity and +infinity With a mean=0.0 and a standard devia- tion=l.0. ******************************************************************** function rnor(j) . patchwork normal, 0