UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Nonradial oscillations in Spica Fraser, Geoffrey Alan 1985

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1985_A6_7 F72.pdf [ 4.86MB ]
Metadata
JSON: 831-1.0085599.json
JSON-LD: 831-1.0085599-ld.json
RDF/XML (Pretty): 831-1.0085599-rdf.xml
RDF/JSON: 831-1.0085599-rdf.json
Turtle: 831-1.0085599-turtle.txt
N-Triples: 831-1.0085599-rdf-ntriples.txt
Original Record: 831-1.0085599-source.json
Full Text
831-1.0085599-fulltext.txt
Citation
831-1.0085599.ris

Full Text

NONRADIAL OSCILLATIONS IN SPICA By GEOFFREY ALAN FRASER B.Sc, The University of V i c t o r i a , 1981 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n THE FACULTY OF GRADUATE STUDIES (Department of Geophysics and Astronomy) We accept this thesis as conforming to the rec|uir#d staadatffl THE UNIVERSITY OF BRITISH COLUMBIA A p r i l 1985 © Geoffrey Alan Fraser, 1985 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. I t i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department Of Geophysics and Astronomy The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date March 15, 1985 -6 (3/81) i i ABSTRACT The absorption l i n e p r o f i l e s of Spica (a V i r g i n i s , HD116658, B1.5IV, m=0.97) show features, at about the 1% l e v e l , moving from the blue wavelengths towards the red wavelengths. A series of spectra were taken, at the 1.22 m telescope at the Dominion Astrophysical Observatory on three nights i n A p r i l , 1982 and two nights i n A p r i l , 198A, to study these moving features. As Spica i s a member of a binary system, the effect of the secondary had to be removed from the observations. This was done by subtracting a template spectrum which had been scaled, broadened and shifted to match the secondary, from each observation. The required s h i f t s were determined using the o r b i t a l elements on blended nights and using the Fahlman-Glaspy small-shifts technique on unblended nights. An average of a l l the spectra was then subtracted from each observation. The resulting series of residuals c l e a r l y show the motion of the features seen i n the l i n e p r o f i l e s . The acceleration of the features was estimated to be between 0.0055 and 0.0068 kms - 2. Assuming the features are due to nonradial o s c i l l a t i o n s , t h i s acceleration corresponds to waves moving slowly, about 5 to 20 kms - 1, i n a prograde d i r e c t i o n . The angular frequency of the o s c i l l a t i o n s , after accounting for the effects of rotation, would be about 3.4X10 - 5 r a d s - 1 . i i i A computer model that produces l i n e p r o f i l e s , under the assumption of a single nonradial o s c i l l a t i o n , was used to produce p r o f i l e s for comparison with observations. Using an 1=8 and m=-8 mode, an i n t r i n s i c frequency of 3.4X10"5 r a d s - 1 and a s t e l l a r rotation rate of 190 kms - 1, the model produced p r o f i l e s similar to those observed. The change i n the model p r o f i l e s with time was also similar to that observed. i v TABLE OF CONTENTS 1.0 INTRODUCTION 1 2.0 THEORETICAL INVESTIGATIONS 5 2.1 Identifying Modes 9 2.2 Rotation Effects 10 2.3 Tidal and Rotational Distortion 12 2.4 Modelling Absorption Line P r o f i l e s 14 2.5 Model Parameters 16 3.0 DATA ACQUISTION AND PROCESSING 24 3.1 Aligning the Spectra 36 3.2 Fahlman-Glsapy Method 39 3.3 Comparison with Calculated Shifts 46 3.4 Secondary Removal 49 3.5 Equivalent Widths 55 3.6 Producing Residuals 63 3.7 Forming the Average 65 3.8 Characteristics of the Residuals 65 3.9 Acceleration of Absorption Features 73 4.0 DISCUSSION 87 References Appendix I 91 93 V LIST OF TABLES TABLE I Observing Parameters 25 TABLE I I Normalisation Polynomials 30 TABLE I I I Orbital Elements Used 47 TABLE IV Comparison between Orbital V e l o c i t i e s Calculated from 48 Alignment Shifts and Those Calculated from Shobbrook's Orbital Elements TABLE V Equivalent Widths of the 4552.654 Si I I I Line 64 TABLE VI Accelerations of Absorption Features i n the Si I I I 4552.6 80 A Line TABLE VII Calculation of the I n t r i n s i c Angular Frequency . 85 v i LIST OF FIGURES Figure 1 A schematic of the Spica binary system. Figure 2 Line p r o f i l e changes as a result of changing X from 6 to 8. Figure 3 The effect of m on l i n e p r o f i l e shape. Figure 4 Line p r o f i l e variations due to an increase i n the maximum ra d i a l v e l o c i t y . Figure 5 The influence of k on the l i n e p r o f i l e . Figure 6 Figure showing 100 pixels of data, centred on the Si I I I 4552.6 A l i n e , as recorded d i r e c t l y after an exposure ( t o t a l spectrum i s 1872 p i x e l s ) . Figure 7 The same data displayed i n Figure 6 after subtraction of the dark exposure. Figure 8 Figure showing the result of applying the correcting polynomials calculated from the f l a t lamp r e s u l t s . Figure 9 A portion of an average lamp. This -lamp i s formed by averaging a number of lamp exposures together. In some cases, i t i s preferable to use a single lamp exposure time, rather than an average lamp, which results i n a signal l e v e l close to the average data l e v e l . Figure 10 The data of Figure 8 after d i v i s i o n by the average lamp of Figure 9. Figure 11 The data of Figure 10 a f t e r f u r t h e r a d d i t i v e l i n e normalization to remove residual video l i n e d i f f e r e n t i a l s . Figure 12 The processed data after f i l t e r i n g to 50% of the Nyquist frequency. Figure 13 A spectrum from A p r i l 21/22, 1982 after processing and r e c t i f i c a t i o n . A t h i r d order polynomial was used to f i t the continuum. The dotted l i n e indicates the normal continuum. Figure 14 Same spectrum as Figure 13 showing the entire observed wavelength range. Figure 15 The 21 observations of the Si I I I t r i p l e t at 4560 A, taken on A p r i l 19/20, 1982. Each spectrum i s labelled with i t ' s midexposure time. The secondary shows up as the weak lines just to the blue of the f i r s t two Si I I I l i n e s . v i i Figure 16 The 10 observations of the Si I I I t r i p l e t at 4560 A, taken 42 on A p r i l 19/20, 1982. Figure 17 The 28 observations of the Si II t r i p l e t taken on A p r i l 21/22, 1982. The secondary appears as the weak lines just to the red of the f i r s t and last Si I I I l i n e s . Figure 18 The 37 observations of the Si I I I t r i p l e t taken on A p r i l 26/27, 1984. The secondary i s completely blended. Figure 19 The 29 spectra of the Si I I I t r i p l e t taken on A p r i l 28/29, 1984. The secondary i s completely blended. Figure 21 The A p r i l 21/22, 1982 data after subtracting the mean of the A p r i l 28/29, 1984 data. The 'S' marks the position of the secondary. Figure 22 The dotted l i n e i s HR6588 and the s o l i d l i n e i s the Si I I I l i n e from Spica. HR6588 has been broadened by convolving i t with a 50 km s - 1 rotation p r o f i l e . Figure 23 The dotted l i n e i s HR6588 after being broadened then s h i f t e d by the amount c a l c u l a t e d from the o r b i t a l elements. 43 44 45 Figure 20 The A p r i l 19/20, 1982 data after subtracting the mean of 51 the A p r i l 28/29, 1984 observations. The 'S' marks the position of the secondary. 52 54 56 Figure 24 The dotted l i n e i s the o r i g i n a l data and the s o l i d l i n e i s 57 the result of subtracting the template spectrum shown i n Figure 23. Figure 25 The A p r i l 19/20, 1982 data set after the secondary has been 58 removed. The obvious variations i n the l i n e p r o f i l e s are due to veloc i t y f i e l d s on Spica. Figure 26 The Si I I I t r i p l e t of Spica after removal of the secondary 59 from data of A p r i l 20/21,1982. Figure 27 The Si I I I t r i p l e t of Spica after removal of secondary from 60 data of A p r i l 21/22, 1982. Figure 28 The Si I I I t r i p l e t of Spica after removal of the secondary 61 from the data of A p r i l 26/27, 1984. Figure 29 The Si I I I t r i p l e t of Spica after removal of the secondary 62 from the data of A p r i l 28/29, 1984. Figure 30 Residuals of the Si I I I 4552.654 lines observed on A p r i l 66 19/20, 1982. Capital l e t t e r s indicate the features whose r a d i a l v e l o c i t y was measured. v i i i Figure 31 Residuals of the Si I I I 4552.654 lines observed on A p r i l 67 20/21, 1982. Figure 32 Residuals of the Si I I I 4552.654 li n e s observed on A p r i l 6 8 21/22, 1982. Figure 33 Residuals of the Si I I I 4552.654 lines observed on A p r i l 26/27, 1984. 69 Figure 34 Residuals of the Si I I I 4552.654 lines observed on A p r i l 70 28/29, 1984. Figure 35 Selected residuals showing dominant pattern seen i n each 72 data set from the f i v e nights of observation and the model. Figure 36 A plot of the measured r a d i a l v e l o c i t y versus time for the 75 four absorption features from the A p r i l 19/20, 1982 data set. Radial velocity was measured r e l a t i v e to the l i n e centre. Figure 37 Radial vel o c i t y versus time for absorption features i n the 76 A p r i l 20/21, 1982 data set. Figure 38 Radial vel o c i t y versus time for absorption features i n the 77 A p r i l 21/22, 1982 data set. Figure 39 Radial vel o c i t y versus time for absorption features i n the 78 A p r i l 26/27, 1984 data set. Figure 40 Radial vel o c i t y versus time for absorption features i n the 79 A p r i l 28/29, 1984 data set. Figure 41 A model time s e r i e s produced using the no n r a d i a l 82 o s c i l l a t i o n model with the following parameters: Veq=190 km s" 1; 1=8; m=-8; k=0.1; V r m a x=8 km s - 1 . Figure 42 The r a d i a l v e l o c i t y versus time for the model residuals. 83 Figure 43 Residuals produced by subtracting the mean of the model 84 time series from each indi v i d u a l model l i n e p r o f i l e . Figure 44 Five model absorption lines versus observed lines from the 88 A p r i l 19/20, 1982 data set. Figure 45 Model residuals plotted versus residuals from the A p r i l 89 19/20, 1982 data set. The phase difference between each model residual corresponds to the time difference between each observed residual. The model phase increments were calculated as (900 sees)* (2.54X10-1* rad s - 1 ) , producing an observed acceleration of the model absorption features of about 0.0055 km s - 1 (see Table VI). i x Acknowledgements I would l i k e to thank my supervisor, Dr. Gordon Walker, for suggesting this topic and for his support and encouragement throughout t h i s investigation. I would also l i k e to thank Stephenson Yang, Chris Millward and Zoran Ninkov, for acquiring several nights of observations including the A p r i l , 1982 data set, Dr. Greg Fahlman and Daniel Thibault, who designed and programmed the algorithm on which the computer model i s based, and Gerry Grieve, who assisted me i n using the VAX computer. My special thanks to John Amor, for his constant aid and his many helpful suggestions, and to Christine Third for her support and help i n assembling this thesis. Financial support was provided by a grant from the National Science and Engineering Research Council. 1 1.0 INTRODUCTION Spica i s an early type star apparently undergoing nonradial o s c i l l a t i o n s (a V i r g i n i s , HD116658, B1.5IV, m=0.97). While other examples of high order o s c i l l a t o r s e x i s t , Zeta Ophiuchi for instance (Vogt and Penrod, 1983), Spica i s unusual i n that i t was a singly periodic 8 Cephei star as recently as 1970. 8 Cephei stars show l i g h t and v e l o c i t y variations of small amplitude and short period. The l i g h t variations are t y p i c a l l y less than 0.1 magnitudes and have periods from three to seven hours. The r a d i a l v e l o c i t y changes are generally less than 50 kms - 1 (Lesh and Aizenman, 1978). Spica was i d e n t i f i e d as a 8 Cephei variable, with a magnitude range of 29.0 m.mag (m.mag = 1 0 - 3 magnitudes) and a 4.17 hour period, i n 1968 (Shobbrook et a l . , 1969). Observations i n 1970 showed the same period but there was a decline In the photometric variations to 14 m.mag(Shobbrook et al,1972). Subsequent observations, i n 1970 and 1971 (Dukes, 1974), showed that four periods, of 4.17, 4.15, 5.83 and 6.63 hours respectively, had appeared i n the vel o c i t y and photometric data. After 1971, these variations declined i n amplitude u n t i l they were undetectable by 1973. Further observations i n the UV (Hutchings and H i l l , 1977) and v i s i b l e wavelengths (Lamb, 1978) revealed no evidence of pulsations. 2 A L P H A VIRGINIS a = 1 . 9 3 * 1 0 7 k m e = 0 . 1 4 P= 4 . 0 1 5 d a y s P a = l 2 4 y e a r s Figure 1 A schematic of the Spica binary system. 3 While the photometric and v e l o c i t y v a r i a t i o n s of Spica e s s e n t i a l l y disappeared by 1972, and have not been observed up to the present (Smith, 1985), another type of v a r i a b i l i t y has been detected. In 1978, Walker et a l . (1982) started a series of observations of Spica. Using a Reticon detector to obtain good time resolution and high signal to noise r a t i o spectra, they demonstrated the presence of absorption l i n e p r o f i l e variations i n the He I 6678 A l i n e . Many B stars show evidence of some v a r i a b i l t y i n their l i n e p r o f i l e shapes. In most cases these variations appear as a changing asymmetry of the l i n e (Smith, 1979). Spica shows much more complex variations than are usually seen. Similar p r o f i l e variations have been observed i n Zeta Ophiuchi (Walker et a l . , 1979). Vogt and Penrod (1983) proposed two models to explain the t r a v e l l i n g features seen i n Zeta Ophiuchi. The f i r s t involves an orb i t i n g structure of occulting material about the star. While this can produce variations similar to those observed i t suffers from at least two problems. In the case of either Spica or Zeta Oph., orbiting material should produce l i g h t variations but none are observed. In addition, the close binary companion of Spica would tend to t i d a l l y disrupt any o r b i t t i n g structure, which would make the observation of persistent features u n l i k e l y . Vogt and Penrod then t r i e d modelling the variations with a high order nonradial o s c i l l a t i o n . This model was quite sucessful i n producing p r o f i l e variations similar to those observed. High order 4 nonradial modes are also unlikely to cause photometric, or equivalent width variations. This i s p a r t i c u l a r l y important given the lack of photometric or equivalent width v a r i a b l i t y i n Spica. The o r i g i n of nonradial o s c i l l a t i o n s i n Spica, or any B star, i s s t i l l uncertain. Many attempts have been made to l i n k the i n s t a b i l i t y of 8 Cephei stars to their position near the S-bend of their evolutionary tracks. Unfortunately, no convincing case can be made for a driving mechanism linked to the evolutionary stage (Unno et a l . , 1979). Driving mechanisms linked to some special characteristic of a star, such as a rapidly spinning core, have also been proposed (Osaki, 1975). The a p p l i c a b i l i t y of these models to Spica remains unclear. Despite questions concerning their o r i g i n , nonradial o s c i l l a t i o n s do produce many of the observed ch a r a t e r i s t i c s of p r o f i l e variables l i k e Spica. I t i s my intention to investigate this behaviour more f u l l y by comparing model l i n e p r o f i l e s , produced under the assumption of the presence of nonradial o s c i l l a t i o n s , with observed p r o f i l e s . I t i s also possible, from comparison with the observations, to gain some idea of the period, and hence the type of nonradial o s c i l l a t i o n that i s occuring. 5 2.0 THEORETICAL CONSIDERATIONS To model the effect of nonradial o s c i l l a t i o n s , i t i s necessary to determine the velo c i t y f i e l d s produced at the star's surface by these o s c i l l a t i o n s . Once these vel o c i t y f i e l d s are known one can generate model spectral absorption lines for comparison with observations. S t e l l a r o s c i l l a t i o n s can be described by a set of nine nonlinear, p a r t i a l d i f f e r e n t i a l equations. These equations describe mass, energy and momentum conservation, f l u x transfer, and the gravitational f i e l d , a l l i n three dimensions. Even i n the r e l a t i v e l y simple case of a spherically symmetric, nonrotating star, where this system of equations can be reduced i n complexity, solutions for a r e a l i s t i c s t e l l a r model are d i f f i c u l t to compute. As a resu l t the further simplifying assumption of small o s c i l l a t i o n s i s generally made. When the o s c i l l a t i o n s are small enough to be considered as pertubations of the undisturbed star, one can substitute variables of the form x + 6x into the system of equations outlined above. The variable x can be any physical quantity, such as position, v e l o c i t y , temperature, etc., and 6x i s the difference from the unperturbed state. Neglecting a l l powers above the f i r s t and products of the variations, results i n a system of line a r d i f f e r e n t i a l equations describing the behaviour of the pertubations. Cox(1980) provides a very detailed look at the entire l i n e a r i z a t i o n procedure. 6 In the l i n e a r , n o n r a d i a l - o s c i l l a t i o n regime, one generally assumes that any scalar physical variable i s the product of a function of r a d i a l d i s t a n c e , a spherical harmonic Y™(9,<|>) and a time dependence e " i a t (Cox, 1980). In turn, Y™(9,<|>) = P™(cos 9 ) e i m * where P°(cos 9) i s a Legendre polynomial of order A and harmonic index m. Index m can take the values -A, -1+1, ... 0, ... 1-1, 1. A small r a d i a l displacement at the surface of a star would therefore be described by: C r= I 4*k(r) ^ ( 9 , * ) e - i a t [2.1] k,l,m ' where k indicates the number of nodes s a t i s f y i n g , 4^k ( ri ) = 0 ( i = i> 2 k ) [ z- z ] In the linear regime, the various o s c i l l a t i o n modes do not interact so that we can drop the summation and consider each mode separately. For a general displacement at the surface we have the r e l a t i o n (Unno et a l . , 1979): - ( y r ) > • V r > i d h H ) ^ <M>e- 1 0 t [2.3] In practise this gives us three displacement equations: 7 x r = -A P™(cos9)sin(m(t)-0t) [2.4] x = -k A 5 P l ( c 0 s 9 ) s i n ( m ( ^ a t ) [2.5] 9 r "o9 x, = —k A m P m (cos9) cos(m<|>-at) [2.6] $ r A sin9 and after taking the time derivative three v e l o c i t y equations: v r = a A r P° cos(m<t)-at) [2.7] v Q = a k A oP m cos(m(|)-0t) [2.8] 9 r 1 59 v = -a k A m P m sin(m<t>-0t) [2.9] (|) r X. sin 9 where 9 i s the polar, and (> the azimuthal angle. The appearance of the star's surface i s governed by the spherical harmonic P™. Moving around the equator of the star, one would encounter m wavecrests, each separated by 2it/ \n \ radians. As each wavecrest represents a particular phase of the disturbance, a crest w i l l propagate with a veloc i t y equal to the phase velo c i t y of the disturbance. Given a time dependance of e *"0t one gets a phase velo c i t y V, given by: 8 V a [2.10] m Thus, for m negative, a wave propagates i n the direction of decreasing azimuth. The angular frequency, a, that characterises a nonradial o s c i l l a t i o n i s set by the restoring force producing the o s c i l l a t i o n . Two restoring forces are possible - pressure or gravity. Pressure variations are characterised by the Lamb frequency, given by: where P = Local pressure p = Density r = Distance from the star's centre Adiabatic constant A disturbance with gravity as the restoring force has a frequency N, given by: A(A+1) r xp [2.11] 1 N = (-A g ) 2 [2.12] where 1 dp _ 1 dP [2.13] p dr I\P dr 9 N i s the Brunt-Uaisala frequency. The factor A indicates s t a b i l i t y (A<0) or i n s t a b i l i t y (A>0) against convection (Cox,1980). I f a 2 of a wave i s greater than N 2 or S 2 then the disturbance i s a pressure or P-wave. For o 2<N 2,S 2 the o s c i l l a t i o n i s a gravity or G-wave. For frequencies between N 2 and S 2 n e i t h e r type of wave can propagate. A region where cr2 for a particular disturbance f a l l s between N 2 and S 2 i s c a l l e d an evanescent region f o r that frequency (Osaki, 1975). G- or P-waves are further characterised by the number of nodes a wave has between the star's centre and surface. In general, the higher the mode the longer the period of a G-wave, and the shorter the period of a P-wave. Thus very long period o s c i l l a t i o n s , of order eight hours or longer, are l i k e l y to be G modes (Osaki, 1975). 2.1 Identifying Modes A number of methods exist to determine whether a star i s o s c i l l a t i n g i n a r a d i a l or nonradial mode. These generally involve comparing observational results with theoretical models. One good method i s to compare observational Q values with theoretical values. Q i s defined by: 10 p I Q = P ( — ) 2 [2.14] where P i s the period and p i s the mean density. For B Cephei stars, theoretical Q values give 0.037 days for the fundamental r a d i a l mode, 0.026 - 0.028 days for the f i r s t harmonic, and 0.022 days for the second harmonic r a d i a l mode. Nonradial modes give Q's of 0.028 days for p^ and 0.022 for p2« When Spica was o s c i l l a t i n g with a 4.17 hour period, i t ' s Q value was about 0.025 which i s consistent with second harmonic r a d i a l modes, or the p^ mode (Lesh and Aizenman, 1978). This thesis, however, w i l l i d e n t i f y the l i n e p r o f i l e variations with a period of about f i f t y hours, resulting i n a Q of 0.3 days. This value for Q i s probably consistent only with a very high order G mode. Unfortunately, no theoretical models of high enough order modes exist for comparison with the observations. 2.2 Rotation Effects When a star i s rotating the frequency of an o s c i l l a t i o n i s changed from i t s value i n the nonrotating case. For an observer i n an i n e r t i a l frame,, observing a star rotating with angular vel o c i t y Q, two 11 effects are important, one due to the c o r i o l i s acceleration and the other due to the frame of reference of the observer. The d i r e c t e f f e c t of r o t a t i o n i s due to the c o r i o l i s acceleration. For a wave t r a v e l l i n g i n the d i r c t i o n of rotation, a prograde wave, the c o r i o l i s force acts to reduce the l o c a l gravity. With reduced gravity, the phase velo c i t y of gravity waves becomes even slower. The opposite effect occurs for waves moving i n a retrograde dire c t i o n (Cox, 1984). For an observer stationary r e l a t i v e to the star, the observed phase velo c i t y of a wave, V ^ g i s equal to: V . - V T M_ + Q [2.15] obs INT 1 1 where V^N^, i s the i n t r i n s i c phase v e l o c i t y perturbed by the c o r i o l i s effect outlined above. Remembering that V = cr/m, we have, a . = aT.Trp + mQ [2.16] obs INT 1 1 For a prograde wave, and mQ have the same sign so that one observes a larger phase ve l o c i t y and hence, a larger angular frequency than otherwise. For a retrograde wave, the observed frequency w i l l be reduced by mQ. For fast rotation, or small i n t r i n s i c v e l o c i t y , a retrograde mode may appear to be dragged i n a prograde direc t i o n by the rotation. 12 Combining the two effects outlined above gives the usual expression: o = a + m(l-C, „)Q [2-17] O K. y X, where C^ ^ i s due to the c o r i o l i s a c c e l e r a t i o n and v a r i e s with the o s c i l l a t i o n mode and the star's structure. There i s a further contribution to the observed angular frequency due to the centrifugal acceleration. However, th i s i s a second order effect and w i l l be ignored for this work. The v a l u e of C^ ^ i s f r e q u e n t l y taken as 0.179 f o r 1=2 o s c i l l a t i o n s of (3 Cephei stars, but i s probably much less than this value f o r higher order G modes (Saio, 1981). Smith (1985) uses Cfc =^ 1 = 0.01 for 1=8. 1(1+1) 2.3 Tidal and Rotational D i s t o r t i o n Deviations from sphericity due to t i d a l and rotational effects can be found using the relations (Shobbrook et al.,1969): M2 . R . , Tidal [2.18] A E = M - ' ( — ) 3 R 3u 2 Rotational [2.19] A T I = "GM7 13 where M-^  = Primary mass = 10.9 M @ M 2 = Secondary Mass = 6.8 M Q R = Radius of primary = 8.1 R Q D = Separation between stars = 1.93x101° km w = Angular vel o c i t y of the primary The rotational d i s t o r t i o n , at less than one percent of the radius, was neglected. Tidal effects cause the equatorial major axis to be distorted by: E q u a t o r i a l Major Axis = 2R(1+Aa) [2.20] 3 where Aa=yAe f o r a Roche model. In t h i s case, Aa=0.023 and the resulting bulge increases the radius by 1.3xl0 8 metres. Assuming a four day o r b i t a l period, one point on the star's surface w i l l traverse this distance i n 24 hours requiring a v e l o c i t y of 1.5 kms - 1. Models incorporating a t i d a l d i s t o r t i o n up to 3 km/s were used with very minor changes to the l i n e p r o f i l e s . I t seems l i k e l y that the t i d a l influence of the secondary may affect the mode excited, but doesn't seriously disturb the resulting p r o f i l e s d i r e c t l y . 14 2.4 Modelling Absorption Line P r o f i l e s The model used i n this work i s quite simple i n p r i n c i p l e . A star's surface i s assumed to consist of a number of discrete surface area elements. About 5400 elements were used i n each model l i n e produced in this thesis. Each surface element i s considered i n turn. A displacement, due to the nonradial mode used, i s calculated for each corner of the unit surface. The angle between the normal to the distorted surface and the l i n e of sight to the observer i s then calculated. I f the surface i s v i s i b l e i t contributes to the resulting l i n e p r o f i l e . About 2900 surface elements contributed to each model p r o f i l e . Every surface element i s assumed to produce an i d e n t i c a l i n t r i n s i c l i n e p r o f i l e of the form exp (-( AX./AX ) 21 with AX = X V, /C. V, e v y o ' o o b b i s termed the broadening v e l o c i t y and was set at 8 kms - 1. The broadening vel o c i t y represents the thermal and microturbulent v e l o c i t i e s of the emitting atoms. The i n t r i n s i c p r o f i l e i s a reasonable model for metal l i n e s , such as Si I I I , i n a B-star atmosphere (Stamford and Watson, 1976). The assumption of invariant i n t r i n s i c l i n e p r o f i l e s provides a great s i m p l i f i c a t i o n i n the model, but i t i s not a very r e a l i s t i c approach. Detailed models of the l i n e p r o f i l e s formed i n s t e l l a r atmospheres have demonstrated that a variety of limb variations occur. In l o c a l thermodynamic equilibrium (LTE) models, a l i n e p r o f i l e tends to 15 vanish at the limb. Non-LTE models produce li n e s that i n i t i a l l y diminish away from the centre but which increase again at the limb. Introduction of vel o c i t y gradients i n the atmosphere can produce large asymmetries i n l i n e p r o f i l e s formed at the centre of the s t e l l a r disk, with smaller asymmetry at the limbs (Mihalas, 1979). Stamford and Watson (1980) considered the effect various p r o f i l e changes would have on a modelled l i n e . They concluded that the i n t r i n s i c l i n e p r o f i l e technique gave results e s s e n t i a l l y i d e n t i c a l to those produced with models incorporating v e l o c i t y gradients and LTE e f f e c t s . Unfortunately the models used i n the Stamford and Watson study included pulsations only up to Jl=3. I t would be a useful exercise to extend this type of study to higher order modes. For the present, however, i t w i l l be assumed that the i n t r i n s i c l i n e p r o f i l e assumption i s at least approximately v a l i d . After determining i f a surface element i s v i s i b l e , the v e l o c i t y at the centre of the surface area i s calculated. The v e l o c i t y i s a combination of the rotation v e l o c i t y at the point being considered and the v e l o c i t y due to the nonradial o s c i l l a t i o n . The i n t r i n s i c l i n e p r o f i l e can then be doppler shifted by the resulting t o t a l v e l o c i t y . As each area element covers the same physical area of the s t e l l a r surface, the contribution from an element must be scaled. This scaling accounts for differences i n the projected surface area, limb darkening and gravity darkening. 16 The projected surface area varies as the cosine of the angle between the surface normal and the l i n e of sight. Thus an area element with i t ' s normal along the l i n e of sight has a projected area equal to i t ' s t o t a l area. If the normal to the area i s at a large angle to the l i n e of sight, the area would contribute l i t t l e to the t o t a l l i n e p r o f i l e . Limb darkening i s accounted for by scaling the surface area with the r e l a t i o n : S = S ( l - 8(l-cos9)) [2.21] where 8 = limb darkening constant and 9 = the angle between the centre of the disk and the centre of the surface area. Gravity darkening i s accounted for by the r e l a t i o n : S = | S * G/G q I [2.22] where G q i s the gravity at the pole and g i s the l o c a l gravity. 2.5 Model Parameters Input to the model consisted of a series of variables describing the s i z e , oblateness, polar gravity, and rotation rate of the star. Information regarding the mode and order, the maximum ra d i a l v e l o c i t y and the apparent phase ve l o c i t y of the o s c i l l a t i o n was also required. 17 Appendix I contains a l i s t i n g of the program from which the order and format of the input can be determined. Figure 2 to 5 show the effect of the various o s c i l l a t i o n parameters. In general terms, the o s c i l l a t i o n mode, 1, determines how complicated the surface v e l o c i t y d i s t r i b u t i o n can become. Figure 2 shows the result of changing 1 from 6 to 8. The 1=8 p r o f i l e shows four d i s t i n c t features while 1=6 has only three. The o s c i l l a t i o n mode m determines how many wavecrests w i l l appear around the star's equator. For m=l, m crests appear, centred on the equator. I f m i s less than 1 then the star's surface i s divided into zones by A— |m | l a t i t u d i n a l boundaries. Each zone w i l l then have m wavecrests. A l i n e p r o f i l e w i l l have the most pronounced distortions for m=l, as shown i n Figure 3. Changing the maximum r a d i a l v e l o c i t y of an o s c i l l a t i o n has a st r a i g h t f o r w a r d r e s u l t . As V i n c r e a s e s , the l i n e d i s t o r t i o n s ° r max increase. The number of features i n a l i n e does not seem to be changed by changing V . Figure 4 shows the result of changing V from o to +20 kms -1 Nonradial o s c i l l a t i o n s have a r a d i a l , polar and azimuthal component. I t was assumed that the r a d i a l componenet was dominant and therefore k was set equal to 0.1 for the model l i n e s . The factor k determines the fraction of V that can be attained by the polar and r max 18 L INDEX TEST : M=-L : Vr=12 Km/S Figure 2 Line p r o f i l e changes as a result of changing & from 6 to 8. 19 M INDEX TEST : L=8 Vr=12 Km/S Figure 3 The effect of m on l i n e p r o f i l e shape. 20 VELOCITY TEST : L=8 M=-8 \]% OF CONT 4548 4552 4556 Wavelength in angstroms Figure 4 Line p r o f i l e variations due to an increase i n the maximum r a d i a l v e l o c i t y . 21 K TEST : L=8 M=-8 : Vr=10 Km/S Figure 5 The influence of k on the l i n e p r o f i l e . 22 azimuthal components of the o s c i l l a t i o n . Figure 5 shows the effect of increasing k from 0.1 to 1.0. In attempting to match model li n e s to observations the most important parameters appear to be A and V r • The o s c i l l a t i o n mode A determines the complexity of the p r o f i l e , while V affects the size v J r max of the features within the p r o f i l e . The o s c i l l a t i o n order m can also change the appearance of the l i n e features, but tends to make the p r o f i l e s too smooth when m i s not equal to A. Obviously the model has a large number of variables that can affect the resulting l i n e p r o f i l e . I f we allow more than one mode to contribute the number of variables i s further increased. One could wonder, t h e r e f o r e , i f a s o l u t i o n , i n terms of A, m, V , rotation ' ' ' ' ' r max' rate and was unique. Indeed, from experiments with the model i t appears that an individual l i n e p r o f i l e can be matched equally well by a number of different choices of A, m, V ' ' r max When trying to match a time series of observations, however, the choice of parameters becomes more r e s t i c t e d . If one t r i e s to match the l i n e residuals the choice of parameters becomes much narrower. A residual i s formed by subtracting the average Spica spectrum from each Individual observation. Given a higher resolution, and a better idea of the impulse response of the telescope-spectrograph system, i t may be possible to narrow down the choice of parameters even further. The 23 parameters chosen as the best f i t to the data i n the thesis are discussed i n the section on r e s u l t s . 24 3.0 DATA ACQUISITION A l l of the spectra used i n t h i s work were obtained at the Dominion Astrophysical Observatory i n V i c t o r i a , B.C.. The coude spectrograph of the 1.22 metre telescope was used with a 1200 line/mm grating giving a dispersion of 10 A/mm. Table I l i s t s the observing system parameters. A 1-D array of 1872 diodes was used as the detector. Each diode i s back biased to f i v e v o l t s . Light, incident on a diode, generates electron-hole pairs which discharge the bias. A readout consists of rebiasing each diode. The current required for each diode i s converted to a voltage pulse which can be integrated, d i g i t i z e d and recorded (see Walker, 1984). A diode array has a number of advantages over the t r a d i t i o n a l photographic plate. These include s p a t i a l s t a b i l i t y , high quantum eff i c i e n c y and nearly lin e a r response. Also, when cooled to l i q u i d nitrogen temperature, the thermal noise i s minimal. There are three problems associated with this detector that must be corrected during the i n i t i a l processing: 1) Four separate c i r c u i t s , or video l i n e s , are used to read out the array. C i r c u i t one reads diodes one, f i v e , nine, etc., 25 TABLE I OBSERVING PARAMETERS Telescope : 1.22 m Spectrograph: 32 inch camera Grating: 1200 lines/mm Dispersion: -10 A/mm (0.15 A/pixel) Resolution: 0.2 - 0.3 A DETECTOR: 1872 Reticon array 26 c i r c u i t two handles diodes two, s i x , ten and so on. During the readout, each of these c i r c u i t s adds noise to the recorded signal. Each c i r c u i t produces a different amount of noise. Readout noise can be corrected for by the use of a dark exposure. A dark exposure consists of an exposure i n which there i s no l i g h t incident on the array. When a dark exposure i s subtracted from the data, the readout noise i s largely removed. As the detector i s not completely stable with time these dark exposures should be taken as close to an observation as possible, and should be of about the same length. During the Spica observations, dark exposures were taken at the beginning and end of each time series. Figure 6 shows an example of raw data while Figure 7 shows the same data after subtraction of a dark exposure. 2) There are diode to diode differences i n response to a given input s i g n a l . These variations are removed by dividing the data by an exposure made with a continuous c a l i b r a t i o n source. The c a l i b r a t i o n source used for the data i n this thesis consisted of a lamp which illuminated the entrance s l i t of the spectrograph. To illuminate the s l i t evenly, the lamp l i g h t i s f i r s t passed through a d i f f u s e r . 3) After subtraction of the dark exposure there are usually small systematic differences between the four video l i n e s . 27 4831.429 4 5 2 3 . 9 6 4 CO C J Q d 4 2 1 6 . 4 9 6 3 9 0 9 . 0 3 3 3601.56( '45 270 (1 (III 295 320 PIXEL # 345 Figure 6 Figure showing 100 pixels of data, centred on the Si I I I 4552.6 A l i n e , as recorded d i r e c t l y after an exposure ( t o t a l spectrum i s 1872 p i x e l s ) . 28 3659.463 PIXEL H Figure 7 The same data displayed i n Figure 6 after subtraction of the dark exposure. 29 These differences show up i n the jagged appearance of the continuum i n Figure 7. To correct for this problem, a series of exposures, of increasing length were taken with a f l a t lamp. This f l a t lamp consists of an ordinary lamp f i l t e r e d to give a nearly constant intensity across the array. Since a l l of the diodes should record an i d e n t i c a l exposure, i t i s possible to calculate a normalisation factor for each video l i n e that brings a l l four video l i n e s into agreement. A polynomial can be f i t t e d to the resulting series of intensity-normalisation pairs for each l i n e . Any exposure, within the intensity range defined by the lamp series, can then be corrected by these polynomials. One series of f i l t e r e d lamps was taken i n 1982, on the night of A p r i l 21/22, and another was done i n May of 1984. The resulting polynomial coe f f i c i e n t s are l i s t e d i n Table I I . A f i n a l additive video l i n e normalisation i s performed on the data to reduce any residual differences. The data was then f i l t e r e d to 50% of the Nyquist frequency. Figures 6 to 12 show the results of each step of the processing flow. A l l data reduction was carried out using the RETICENT processing package at UBC (Yang, 1982). Following the i n i t i a l processing steps, each spectrum was r e c t i f i e d . This was done by dividing the spectrum by a polynomial -30 TABLE I I NORMALIZATION POLYNOMIALS Data Set Line No. A Coefficients B C A p r i l 1982 1 1.0040195 2.5759X10-8 2 0.98561768 2.20277X10-8 3 1.0220751 3.9429X10-8 4 0.98815544 -1.39704X10-7 A p r i l 1984 1 1.0015667 2 0.98834617 -1.21824X10"6 2.176682X10-10 3 1.0212082 4 0.9888845 1.375076X10"6 -2.6926285X10"10 The coefficients A, B, and C are defined as follows: Y = A + Bx + Cx 2 + ... 31 Figure 8 Figure showing the result of applying the correcting polynomials calculated from the f l a t lamp r e s u l t s . 32 1 .315 1 243 1 1 1 1 1 1 245 270 295 320 345 PIXEL # Figure 9 A portion of an average lamp. This lamp i s formed by averaging a number of lamp exposures together. In some cases, i t i s preferable to use a single lamp exposure time, rather than an average lamp, which results i n a signal l e v e l close to the average data l e v e l . 33 Figure 10 The data of Figure 8 after d i v i s i o n by the average lamp of Figure 9. 34 2846.753 PIXEL n Figure 11 The data of Figure 10 a f t e r f u r t h e r a d d i t i v e l i n e normalization to remove residual video l i n e d i f f e r e n t i a l s . 35 2846.406 PIXEL tf Figure 12 The processed data after f i l t e r i n g to 50% of the Nyquist frequency. 36 fourth order for 1982 data and third order for 1984 data, f i t t e d to selected continuum points. Figures 13 and 14 show the results for the 4552 A Si I I I l i n e and the entire spectrum, respectively. 3.1 Aligning the Spectra Spica i s the primary of a binary system. As a r e s u l t , each spectrum i n a time series i s shifted i n wavelength r e l a t i v e to the previous spectrum. Two approaches to correct for this s h i f t were considered. If the elements of the binary orbi t are known, the expected r a d i a l v e l o c i t y of Spica can be calculated for any given time. Corresponding wavelength changes can be found and the data shifted accordingly. This approach requires that good o r b i t a l elements be available. Two sets of elements were considered, those due to Shobbrook et a l . (1972) and a set from Moyles (1982). There are a number of d i f f i c u l t i e s i n determining an orbit for Spica. The primaries' spectral lines are broad (about Vsin i of 160 km s - 1 ) and of variable shape. Spectral l i n e s due to the secondary are blended with the primary for at least half the four day o r b i t a l period. There i s also some evidence that the apsidal period may be changing with time (Moyles, 1982 ). 37 1 .002 PIXEL « Figure 13 A spectrum from A p r i l 21/22, 1982 after processing and r e c t i f i c a t i o n . A third order polynomial was used to f i t the continuum. The dotted l i n e indicates the normal continuum. Figure 14 Same spectrum as Figure 13 showing the entire observed wavelength range. 39 Given these kinds of problems i t was unclear which, i f any, o r b i t a l parameters would give the best r e s u l t s . I therefore decided to calculate the required s h i f t s d i r e c t l y using a method developed by Fahlman and Glaspy (1973). 3.2 Fahlman-Glaspy Method This method assumed that two observations of a l i n e , a(X) and b(X), are described by the same i n t r i n s i c p r o f i l e p(X), with b(X) s h i f t e d an amount AX . I f one applies a small known s h i f t AX to b(X), o s a difference function can be constructed: d(X, AXQ, AXg) = a(X) - b(X-AX g) [3.1] or d(X, AX , AX ) = p(X) - p(X-AX +AX ) [3.2] v ' o s s o Expanding the la s t term i n equation 3.2 as a Taylor series , keeping only the lowest order terms, gives: p(X-AXg+AXo) = p(X) - p'(X)AX g + p«(X)AXo [3.3] Substituting into equation 3.2 yields: d(X, AX , AX ) = (AX - AX ) p'(X) v ' s o s O [3.4] 40 To lessen the method's s e n s i t i v i t y to scaling differences between a(X) and b(X), and to noise, one defines a function: y(AX ) = £ d 2(X, AX AX ) [3.5] S X o s where one sums over the entire l i n e p r o f i l e . To f i r s t order, this gives us: y(AX ) = (AX 2 - 2AXJAX + AX 2) S [p'(X)] 2 [3.6] 3 s d d o o . A. The f u n c t i o n y(AX ) i s a parabola with a minimum when AX =AX . I f the 3 s s o li n e p r o f i l e s were i d e n t i c a l , the minimum would be zero. There are several problems i n applying this technique to real observations of Spica. F i r s t of a l l , the i n t r i n s i c l i n e p r o f i l e s are not always symmetric, and they vary with time. Further, when the secondary spectrum i s blended with the primary, the resulting p r o f i l e i s weighted to the wing containing the secondary. Ignoring these problems for the moment, the spectra were aligned using the s h i f t s calculated with the Fahlman-Glaspy technique (as implemented i n the RETICENT package). Figures 15 to 19 show the results for the fi v e time series obtained. To the eye, the alignments look good, although experience indicates that misalignments of one pi x e l or less are d i f f i c u l t to pick out. The s h i f t s applied varied between a A l 4 5 4 0 4 5 4 8 4 5 5 6 4 5 6 4 4 5 7 2 4580 Uavelength 5X OF COST Figure 15 The 21 observations of the Si I I I t r i p l e t at 4560 A, taken on A p r i l 19/20, 1982. Each spectrum i s labelled with i t ' s midexposure time. The secondary shows up as the weak l i n e s just to the blue of the f i r s t two Si I I I l i n e s . 42 Figure 16 The 10 observations of the S i I I I t r i p l e t at 4560 A, taken on A p r i l 19/20, 1982. 43 msxr SI OF COMT 4549 4558 4564 4372 4580 U m Length Figure 17 The 28 observations of the SI II t r i p l e t taken on A p r i l 21/22, 1982. The secondary appears as the weak lin e s just to the red of the f i r s t and last Si III l i n e s . 44 187140 51 OF CONT Figure 18 The 37 observations of the Si I I I t r i p l e t taken on A p r i l 26/27, 1984. The secondary i s completely blended. 45 194353 035752 5* OF CONT 4540 Figure 19 The 29 spectra of the Si III t r i p l e t taken on A p r i l 28/29, 1984. The secondary i s completely blended. 46 high of 4.0 pixels for the last observation of A p r i l 26/27 to a low of about 0.1 p i x e l s . 3.3 Comparison with Calculated S h i f t s Shifts from the Fahlman-Glaspy technique can be converted to v e l o c i t i e s , and compared with v e l o c i t i e s calculated from the o r b i t a l elements (Appendix I contains a l i s t i n g of a program to convert o r b i t a l elements into r a d i a l v e l o c i t i e s ) . As the Si I I I lin e s are quite asymmetric i t i s d i f f i c u l t to determine the centre of a l i n e , and hence the absolute v e l o c i t y . Therefore, I assumed that the ve l o c i t y of the f i r s t spectrum i n a time series equals the calculated v e l o c i t y . Table I I I l i s t s the o r b i t a l elements used, the range i n residuals for each set of observations appears i n Table IV. The residuals were calculated by subtracting v e l o c i t i e s from the s h i f t s from the v e l o c i t i e s calculated using o r b i t a l elements. Several results are apparent. On the two unblended nights, A p r i l 19/20 and 21/22,1982, the mean of the residuals i s about 1.5 km/s. At a dispersion of about 10 km/pixel (.15 A/pixel), t h i s implies a mean difference of about 0.2 pi x e l s . On nights with a blended spectra the mean residuals vary between 3 and 8 km/s, a difference of 0.3 to 0.8 pi x e l s . 47 TABLE I I I ORBITAL ELEMENTS USED (km s" 1) 0 ±2 k l (km s - 1 ) 124 ±4 (km s - 1 ) 197 +8 e 0.14 ±0.03 P (days) 4.01454 +0.00003 O) (°) 142 ±8 T (j.d.) 2440284.78 ±0.08 Q (years) 128 ±12 System velo c i t y Amplitude of the o r b i t a l v a r i a t i o n of the primary Amplitude of the variation of the secondary Time of periastran passage Longitude of periastran at T Orbital e c c e n t r i c i t y Period of o r b i t a l motion Period of rotation of the l i n e of apsides Elements due to Shobbrook et a l . , 1972 48 TABLE IV Comparison between o r b i t a l v e l o c i t i e s calculated from alignment s h i f t s and those from Shobbrooks' o r b i t a l elements Dates Mean Difference (km s" 1) Range of Differences (km s - 1 ) Blended Notes A p r i l 19/20 1982 1.44 -0.05 •* -2.9 U 1 A p r i l 20/21 1982 -3.11 -0.9 * -6.2 B 2 A p r i l 21/22 1982 -1.64 -0.5 * -4.2 U 3 A p r i l 26/27 1984 4.8 -0.07 + -10.9 B 4 A p r i l 28/29 1984 -7.7 -2.7 -»- -9.8 B 5 U - Secondary Blended B - Secondary Unblended 1 Secondary i n the blue wing of the primary and moving to the blue 2 Secondary moving from the blue into the red 3 Secondary in the red wing and moving further into the red 4 Secondary in the l i n e centre at the start of the night, moving toward the red 5 Secondary moving from blue to red across the l i n e centre 49 The residuals are also consistently less than or greater than zero, depending on the location of the secondary. On A p r i l 26/27 1984, the secondary i s moving from the l i n e centre towards the red wing. This appears to weight the effective l i n e centre towards the red. As the primary i s moving blueward this results i n s h i f t s consistently smaller than they should be and hence positive residuals. On A p r i l 28/29, the secondary moves from the blue wing to the red. Thus the i n i t i a l spectrum has i t s centre weighted toward the blue with subsequent spectra weighted more toward the red. This results i n larger than needed s h i f t s and a negative residual. Similar problems are apparent i n the other nights. From these results i t appears that the o r b i t a l elements technique gives s h i f t s good to at least 1.5 kms - 1. As the Fahlman-Glaspy technique has obviously been affected by the secondary on a l l f i v e nights, the calculated s h i f t s may be better than 1 kms - 1. The Fahlman-Glaspy method could be improved on unblended nights by more careful selection of pixel ranges. I t appears, however, that the o r b i t a l elements of Shobbrook et a l . (1972) give very good re s u l t s . 3.4 Secondary Removal As Spica i s the primary of a close binary system, i t ' s spectrum contains a contribution from the secondary. The secondary contribution 50 should be removed to ensure that spectral l i n e variations are s t r i c t l y due to events occuring on Spica. I chose to remove the secondary spectrum by scalin g , s h i f t i n g and broadening a template spectrum to f i t the secondary. Spica's companion i s of the spectral type B3V. I decided to use the spectrum of HR6588 of type B3IV as the template. The f i r s t step i n the procedure was to ali g n the mean of the A p r i l 28/29, 1984 data with the primary spectrum of the A p r i l 21/22 and A p r i l 19/20, 1982 data sets. Required s h i f t s were calculated using the Fahlman-Glaspy method. The A p r i l 28/29, 1984 mean was used since the secondary star's r a d i a l v e l o c i t y d i f f e r s by less than 45 kms - 1 from that of the primary. As a r e s u l t , the secondary should not introduce any spurious asymmetry i n the mean of this night's observations. The mean was then scaled to match each spectrum i n the data set. Subtracting the mean from a l l the observations taken on these two nights l e f t the secondary plus the l i n e p r o f i l e variations. Figures 20 and 21 show the results of th i s procedure. As can be seen from the Figures, the secondary i s almost the same size as the absorption features i n the l i n e , about 1% of the continuum. HR6588 i s a very narrow lined star, V s i n i i s about 13 km s - 1 . This c l e a r l y i s too slow for the secondary l i n e s . To match HR6588 to the secondary, I convolved the observed HR6588 spectrum with a rotation p r o f i l e G(A\) given by (Gray, 1976): 51 11 OF CONT 4 5 4 4 4 5 4 8 4552 4 5 5 8 Uavelenqth In a n g s t r o i s Figure 20 The A p r i l 19/20, 1982 data after subtracting the mean of the A p r i l 28/29, 1984 observations. The 'S' marks the position of the secondary. 52 L8C65B Leen H826B4 0.83718 0J4771 0.92318 133717 A A • V » .« rx OF COHT Figure 21 The A p r i l 21/22, 1982 data after subtracting the mean of the A p r i l 28/29, 1984 data. The 'S' marks the position of the secondary. 53 G(A\) = ^ ( l - C ^ - ) 2)2 + C 2 ( l - ( - ^ - ) 2 ) [3.7] where A\]L=( \ V s i n i ) / c , C1=2 (l-e)/-n;AX 1(l-e/3) and C2=(n;e/2)/'n;A\L(l-e/3) and we assume a limb darkening law of the form: 1 = 1 ° f ( l - e ) + ecosG C C [3.8] where I = the intensity of the continuum and e=0.6. c J Figure 22 shows HR6588 broadened with a 50 kms - 1 p r o f i l e compared with an observed Si I I I l i n e . A broadening velo c i t y of 50 kms - 1 was adopted although v e l o c i t i e s as large as 60 kms - 1 matched some of the l i n e s . After broadening the template, the s h i f t s required to align the Si I I I l i n e s of the template spectrum with the secondary Si I I I l i n e s were calculated. For the A p r i l 19/20 and 21/22, 1982 data, where the secondary was c l e a r l y v i s i b l e , the Fahlman-Glaspy technique was used to a l i g n the two spectra. On nights with the secondary and primary well blended, the position of the secondary was calculated using o r b i t a l elements. A spline f i t was then made between the o r i g i n a l broadened template l i n e positions and the shifted positions. This allowed for the fact that the required s h i f t s varied from one spectral l i n e to the next. The template l i n e s were then scaled, to match peak to peak with the 54 Figure 22 The dotted l i n e i s HR6588 and the s o l i d l i n e i s the SI I I I l i n e from Spica. HR6588 has been broadened by convolving i t with a 50 km s - 1 rotation p r o f i l e . 55 secondary, and shifted to the required positions. Shifts were done using an interpolation routine called CWAVS i n the RETICENT package. Figure 23 shows HR6588 broadened, shifted and scaled versus the f i r s t spectrum of the A p r i l 19/20, 1982 data set. Figure 24 shows the A p r i l 19/20 spectrum after HR6588 has been subtracted, e f f e c t i v e l y removing the secondary. The broadened, scaled and shifted template was subtracted from the data. On nights with the secondary well blended, an average scaling factor from the unblended nights was used. Figures 25 to 29 show the fiv e data sets with the secondary removed. 3.5 Equivalent Widths With the secondary removed from the observed spectra, any changes in the absorption lines should be due solely to the events i n the l i n e forming region of Spica. I f the observed p r o f i l e variations, evident i n Figures 25 to 29 are due to vel o c i t y f i e l d s i n the l i n e forming region, the equivalent widths of the l i n e s should be nearly constant. This w i l l hold true as long as the temperature i n the l i n e forming regions i s not greatly altered by the disturbing o s c i l l a t i o n s . An automated procedure, using the image processing system at UBC, was used to measure the 125 spectra comprising the fi v e data sets. The f i r s t spectrum from a data set was displayed, and the l i m i t s of the l i n e PIXEL H Fieure 23 The dotted l i n e i s HR6588 after being broadened then Figure T h e ^ d ^ ^ c a l c u l a t e d f r o m t h e o r b i t a l elements. 58 Figure 25 The A p r i l 19/20, 1982 data set after the secondary has been removed. The obvious variations in the li n e p r o f i l e s are due to veloc i t y f i e l d s on Spica. 59 Figure 26 The Si I I I t r i p l e t of Spica after removal of the secondary from data of A p r i l 20/21,1982. 60 Figure 27 The Si I I I t r i p l e t of Spica after removal of secondary from data of A p r i l 21/22, 1982. 61 5 1 OF COST 4 5 4 0 4 5 4 8 4 5 5 6 4 5 6 4 4 5 7 2 U j v v l e n q t t i n a n g s t r o a s 4580 Figure 28 The SI I I I t r i p l e t of Spica after removal of the secondary from the data of A p r i l 26/27, 1984. 62 S X OF C Q N T 4 5 4 0 4 5 4 8 4 5 5 6 4 5 6 4 4 5 7 2 4 5 8 0 U a v e L e n q t h in d n g s t r a » s Figure 29 The Si I I I t r i p l e t of Spica after removal of the secondary from the data of A p r i l 28/29, 1984. 63 to be measured picked by eye. Only the 4552 A Si I I I l i n e was measured as i t was least affected by blending. As the variations i n equivalent width, rather than i t s absolute value, are of i n t e r e s t , I decided to include only the core of the l i n e within the l i m i t s chosen. This choice excludes the wings that seem to l i e to either side of the 4552 A l i n e . The area between the chosen l i m i t s was then integrated using a Simpson's Rule algorithm. Table V gives the result for each of the data sets. In general, the equivalent width varied by less than 3% from one spectrum to another. I t seems l i k e l y that more accurate measurements would reduce the variation further. 3.6 Producing Residuals A look at any of Figures 25 to 29 shows that the absorption l i n e p r o f i l e s change shape s i g n i f i c a n t l y over one exposure time, about 900 seconds. A vel o c i t y f i e l d i n the l i n e forming region could redistribute the flux within a l i n e producing variable l i n e p r o f i l e s . An average of such disturbed p r o f i l e s should produce a spectrum corresponding to a v e l o c i t y f i e l d averaged to zero. I f one then subtracts this average spectrum from an individual spectrum, the residual should represent the effect of the v e l o c i t y f i e l d present at the observation time. This procedure was carried out on the f i v e data sets presented i n this thesis. 64 TABLE V Equivalent Widths of the 4552.654 S i I I I l i n e Date Mean Ew (A) Mean Difference % No of Spectra A p r i l 19/20 1982 0.244 0.0052 2.0 21 A p r i l 20/21 1982 0.239 0.0012 0.5 10 A p r i l 21/22 1982 0.235 0.0041 2.0 28 A p r i l 26/27 1984 0.239 0.006 2.5 37 A p r i l 28/29 1984 0.225 0.0026 1.0 29 65 3.7 Forming the Average Ideally an average would be formed from spectra covering many cycles of the disturbing o s c i l l a t i o n s . As more observations are obtained this ideal average w i l l be better approximated. For now, however, only the spectra from the f i v e data sets, plus a few from one night i n May of 1984, were available. As each night's observations have a s l i g h t l y different dispersion r e l a t i o n , the average of each data set was shifted to a l i g n i t with the average of the A p r i l 26/27, 1984 data. Shifts were found using the Fahlman-Glaspy method. Shifted averages were then combined to form a grand-average spectrum. This grand-average was then aligned with each data set and subtracted from the observed spectra. Figures 30 to 34 show the residuals resulting from this procedure. 3.8 Characteristics of the Residuals As the residuals were formed by subtracting the mean from the data, a trough i n the residuals corresponds to a point where the observed spectrum had greater absorption than the mean. These troughs are referred to as absorption features. Several characteristics of the absorption features can be seen by inspecting Figures 30 to 34: 1) Features of approximately constant size and shape persist from one observation to another. 66 1 « OF C O N T 4 5 4 8 4 5 5 2 4 5 5 6 U a v e L e n g t h i n a n g s t r o i s Figure 30 Residuals of the Si III 4552.654 lines observed on A p r i l 19/20, 1982. Capital l e t t e r s indicate the features whose r a d i a l v e l o c i t y was measured. 67 A B C A D 0.86277 0.87330 0.88383 0.89436 0.90489 0.91543 0.92597 0.93650 0.94877 0.96277 1% OF CONT 4 5 4 8 4 5 5 2 4 5 5 6 U a v e i e n g t h i n a n g s t r o m s Figure 31 Residuals of the Si I I I 4552.654 lines observed on A p r i l 20/21, 1982. 68 0 . 8 8 3 8 3 0.30038 0.33717 0.35117 0.36518 0.37918 0 . 9 3 3 1 3 4 5 4 8 4 5 5 2 4 5 5 6 U a v e l e n g t h In a n g s t r o a s IX OF CONT Figure 32 Residuals of the Si III 4552.654 l i n e s observed on A p r i l 21/22, 1982. 69 0.67140 0 . 8 2 9 1 8 0 . 8 3 6 2 4 0 . 8 4 3 3 0 0 . 3 4 9 6 6 0 . 8 5 7 * 2 0 . 3 6 4 1 8 0.87154 0 . 8 7 8 6 0 0 . 8 8 S 6 6 0.89271 0.83631 0.36411 | |11 OF CONT 0.37188 I 4 5 4 8 4 5 5 2 - 6 5 6 U a » e l s f i q t h i n a n q s t r a M Figure 33 Residuals of the Si I I I 4552.654 l i n e s observed on A p r i l 26/27, 1984. 70 0:709 75 0.73967 0.92952 0.94353 0.95 752 4548 4552 4556 U a v e i e n g t h In a n g s t r o » s Figure 34 Residuals of the Si III 4552.654 l i n e s observed on A p r i l 28/29, 1984. 71 2) These persistent features, marked by c a p i t a l l e t t e r s i n the figures, tend to move from the blue wavelengths towards the red. 3) On three nights, A p r i l 19/20 & 20/21, 1982 and A p r i l 28/29, 1984, three persistent features are evident. Only two strong features are seen i n the A p r i l 21/22, 1982 and A p r i l 27/27, 1984 data. Figure 35 shows six residuals, one from each night's observations plus one from a model data set. I consider these particular examples to be representative of the dominant residual pattern of each data set. The persistence of particular absorption features implies an equally persistent o r i g i n . A constant phase point, such as a wavecrest, on a nonradial o s c i l l a t i o n wave could provide the required s t a b i l i t y . A feature's position within the l i n e p r o f i l e would then be set by the r a d i a l v e l o c i t y of the feature's o r i g i n point. To put the motions of the various features on a more quantitative basis, their positions were measured r e l a t i v e to the centre of the l i n e . The l i n e centre was determined by finding the minimum of a parabola f i t t e d to the f i r s t observation of each data set. As each observation i s defined by only three to four p i x e l s , they sometimes assume very asymmetric shapes. Given the few points involved, i t i s d i f f i c u l t to determine a r e l i a b l e centre, for an absorption feature, using a curve 72 10.5% OF CONT 4550 4555 PIXEL Figure 35 Selected residuals showing dominant pattern seen i n each data set from the f i v e nights of observation and the model. 73 f i t t i n g routine. As a r e s u l t , I decided to choose the centre of the features by eye. Measurements were performed on the image processing system at UBC. A residual would be displayed and a cursor moved to the estimated centre of a feature. The position measured could then be compared with the l i n e centre determined e a r l i e r . The miminmum increment of the cursor position was about 0.3 p i x e l s . Errors i n choosing the centre of a feature probably amount to at least this value. Once the position, r e l a t i v e to the l i n e centre, of an absorption feature i s known, i t s r a d i a l v e l o c i t y can be calculated via the usual doppler equation. 3.9 Acceleration of Absorption Features The linear v e l o c i t y of any point on the surface of a rotating star i s given by: - > - » • - > V = Q X R [3.9] where Q = a n g u l a r v e l o c i t y and R = radius vector to any point i n question. The r a d i a l v e l o c i t y corresponds to the projection of this v e l o c i t y along the l i n e of sight. Assuming the l i n e of sight coincides with the z-axis gives: v = x Q s i n i z [3.10] 74 Here x varies fom zero to one, zero being at the star's centre. The corresponding wavelength s h i f t i s given by (Gray, 1976): V X Q s i n i AX = -& = f W [3.11] c v c ' J Thus a l l points on a star's surface having the same x coordinate show the same wavelength displacement and r a d i a l v e l o c i t y . If an absorption feature corresponds to a p a r t i c i l a r o s c i l l a t i o n wavecrest, then a change i n the feature's r a d i a l v e l o c i t y corresponds to a change i n the wavecrest's x coordinate. More precisely, Av Ax Z (0 s i n i ) [3.12] At At where Av /At i s termed the acceleration of a feature and Ax/At i s the z phase ve l o c i t y of the wave. Figure 36 to 40 show the measured v e l o c i t i e s versus time for the various features of Figures 30 to 34. Accelerations of the features, determined by linear regression of the time-velocity data, appear i n Table VI. The measured accelerations range from 0.0025 to 0.0068 kms - 2. As the wavcrests move over a spherical surface, their motions along the x-axis w i l l vary, decreasing as the wavecrests approach the limb. One might expect, therefore, that the measured accelerations would be least for features at large r a d i a l v e l o c i t i e s . This i s true i n 75 VELOCITY RELATIVE TO LINE CENTRE 0 4 - 1 9 - 8 2 2 0 . 0 T 44.0 T 6 8 . 0 9 2 . 0 1 1 6 . 0 , 143.0 TIME (SEC) (X10 2 ) 164.0 188.0 212.0 Figure 36 A plot of the measured r a d i a l v e l o c i t y versus time for the four absorption features from the A p r i l 19/20, 1982 data set. Radial velocity was measured r e l a t i v e to the l i n e centre. 76 VELOCITY RELATIVE TO LINE CENTRE 0 4 - 2 0 - 8 2 8 J • • O D . 1 A A V * ^ A C / r \ v j i • I o o o o o x x x X ^ x x X * X A S _ l 1 1 1 1 1 1 i r 10.0 2 4 . 0 38.0 5 2 . 0 66.CI 8 0 . 0 94.0 108.0 122.0 TIME (SEC) IXIO 2 ) Figure 37 Radial v e l o c i t y versus time for absorption features i n the A p r i l 20/21 data set. 77 VELOCITY RELATIVE TO LINE CENTRE 0 4 - 2 1 - 8 2 8 . a-o s 8 8 GO 0 0 ° ^ O B o o ° A A A A A A A A A * * * X * * * * * x x * w X X A X x X X A x x * x X * * X * X X X 1 1 1 1 1 1 1 1 - 4 . 0 2 6 . 0 SB.O 8 6 . 0 118.0 1 4 6 . 0 , 116.0 206.0 2 3 6 . 0 266 TIME (SEC) (X10 2 ) e 38 Radial v e l o c i t y versus time f o r absorption features i n the A p r i l 21/22, 1982 data set. 78 VELOCITY RELATIVE TO LINE CENTRE 0 4 - 2 6 - 8 4 T 1 1 1 1 1 1 r . 0 31 .0 8 6 . 0 101.0 136.0 171.0, 206.0 241.0 276.0 311.0 TIME (SEC) (X102 ) Figure 39 Radial v e l o c i t y versus time for absorption features i n the A p r i l 26/27, 1984 data set. 79 VELOCITY RELATIVE TO LINE CENTRE 0 4 - 2 8 - 8 4 o si. o o _) mo "8. A A A A A A ^ AA A C A A O O O O O O B m n ^  w w o o O y X X X w >^ X X X a x x * * * x * X x x XXX X X x x x x 5? I 206.0 -1 236.0 2 6 . 0 S 6 . 0 —1 66.0 116.0 TIME (SEC) 1 4 6 . 0 , (X10 2 176.0 266.0 Figure 40 Radial v e l o c i t y versus time for absorption features i n the A p r i l 28/29, 1984 data set. 80 TABLE VI Accelerations of Absorption Features In the Si I I I 4552.654 A Line Date Feature No. of Acceleration Radial Velocity Spectra Range (kms - 2) (kms - 1) A p r i l 19/20 1982 A 21 0.0068±0.0002 -138 to -10 B 17 0.006810.0003 -2 to 90 C 7 0.0052±0.0009 6 to 40 D 21 0.0025±0.0003 120 to 174 A p r i l 20/21 1982 A 10 0.0028±0.0006 -120 to -96 B 10 0.0034±0.0004 -50 to -17 C 10 0.0047±0.0005 20 to 56 D 10 0.0037±0.0002 144 to 174 A p r i l 21/22 1982 A 28 0.0055±0.0002 -134 to 0 B 21 0.0050±0.0001 57 to 144 C 15 0.0057±0.0002 -11 to 52 A p r i l 26/27 1984 A 20 0.0046±0.0002 -110 to -54 B 37 0.0052±0.0002 50 to 175 B 2 19 0.0068+0.0001 47 to 140 A p r i l 28/29 1984 A 29 0.0054±0.0002 -102 to -3 B 21 0.0034±0.0002 46 to 100 C 12 0.0020±0.0003 130 to 144 MODEL 1=8, m=-8 A 25 0.00455±0.00007 -41 to--140 B 25 0.00550±0.00005 -40 to 77 C 25 0.0035+0.0001 81 to 160 81 some instances (feature D, A p r i l 19/20, 1982; features A & D, A p r i l 20/21, 1982; feature C, A p r i l 28/29, 1984) but there are also cases of small accelerations at moderate r a d i a l v e l o c i t i e s . Some indication of what i s happening can be found from the results of the model. Figure 41 shows a model time series produced assuming V =8 kms" 1, 1=8, m=-8, V = 190 kms - 1 and a = 3.42X10 - 5 r max eq o r a d s - 1 . Here i s the equatorial ve l o c i t y of rotation and O"q i s the angular frequency of the o s c i l l a t i o n . An equatorial v e l o c i t y of 190 kms - 1 w i t h an i n c l i n a t i o n of 65.9° (Herbison-Evans, 1971) gave the best match of model l i n e widths to observed l i n e widths. This value of V eq w i l l be assumed from now on. Figure 42 shows the velo c i t y versus time of the residual features produced by subtracting the average model l i n e p r o f i l e from each indiv i d u a l model spectrum. Figure 43 shows the model residuals. As can be seen from Table VI, the measured accelerations from the model data vary from 0.0035 to 0.0055 kms - 2. Only the larger value approaches the acceleration one expects to measure with the given O q . I t would appear, therefore, that the larger accelerations give the best estimate of a wavecrest's actual motion over the star. Using equation 3.12 and the physical characteristics from Figure 1, one can calculate the phase velo c i t y of the o s c i l l a t i o n . Table VII l i s t s the results of these calculations. Once the phase ve l o c i t y i s known the apparent angular frequency can be determined using: 82 2% OF CONT 4 5 4 8 4 5 5 2 4 5 5 6 Uavelength i n angstrons Figure 41 A model time s e r i e s produced using the n o n r a d i a l o s c i l l a t i o n model with the following parameters: Veq=190 km s -1. A=8; m= -8; k-0.1; V r m a x=8 km s~K 83 VELOCITY RELATIVE TO LINE CENTER MODEL o 8 J A A A A A A A A A A A A A A A A A A * 5H x o o 11 | Q o o o X X « x x x x w x x X x x x 8. 1 1 1 1 1 1 1 1 - 4 . 0 2 6 . 0 S 6 . 0 8 6 . 0 116.0 1 4 6 . 0 , 176.0 206.0 238.0 TIME (SEC) (X10 2 ) 5.0 Figure 42 The r a d i a l v e l o c i t y versus time for the model residuals. 84 IX OF CONT 4 5 4 8 4 5 5 2 4 5 5 6 W a v e l e n g t h i n angstroms Figure 43 Residuals produced by subtracting the mean of the model time series from each individual model l i n e p r o f i l e . 85 TABLE VII CALCULATION OF THE INTRINSIC ANGULAR FREQUENCY a f (kms - 1) ±0.0002 V wave (kms - 1) ±10% ( r a d s - 1 ) ±10% a a ( r a d s - 1 ) ±10% a o (r a d s - 1 ) ±10% Period (hrs) 0.0055 179 3.17X10-5 2.54X10-1* 3.2X10_lt 54 0.0057 185 3.28X10-5 2.63X10-1* 4.1X10 - 5 42 0.0068 220 3.90X10-5 3.12X10~'t 9.0X10 - 5 19 190 kms - 1 8 ' 1 R0 65.9° 8 NOTE: 1) Azimuth i s defined as increasing i n the counterclockwise d i r e c t i o n . 2) Spica rotates clockwise so that a prograde wave has a negative phase v e l o c i t y . 3) Since Spica's rotation i s also defined as negative i n this system one takes the difference between a and mQ(l-c) to find cL a . o Using V eq R i m 86 a V = — [3.13] (j) m where V^= phase ve l o c i t y i n rads--'-, o*a= apparent angular frequency, and m = o s c i l l a t i o n mode. Recalling equation 2.17 we have, a = a + mQ(l-c) c~0.1 [3.14] C Si The r e s u l t i n g i n t r i n s i c angular frequency of the o s c i l l a t i o n , a^, i s given in Table VII. 87 4.0 DISCUSSION Since the rotation rate of Spica i s about 190 kms - 1, giving a projected (Vsini) rate of 173 kms - 1, the o s c i l l a t i o n waves are prograde, but not by much. With a projected vel o c i t y between 179 and 220 kms - 1, the waves outpace the rotation by at most 47 kms - 1. Further observations, preferably at higher resolution, may give more consistent values for the acceleration of the features. From Table VII, i t appears that the o s c i l l a t i o n s have a period between about 20 and 50 hours. Myron Smith (1985) has determined a period of 50 hours for an A=16 mode he detected i n Spica's p r o f i l e v ariations. He also finds evidence for an Z=8 mode with a frequency that i s half that for the A=16 mode. Using an i n t r i n s i c period of 54 hours (o"o=3.2X10-5 r a d s - 1 ) and a maximum r a d i a l v e l o c i t y of 8 kms - 1, the nonradial model produced the time series of Figure 41. Some of these model li n e s are compared with observed data i n Figure 44. While the model lines do not reproduce every bump i n the p r o f i l e , they do match the major features. Perhaps more convincing i s the comparison between the corresponding model and observed residuals, as shown i n Figure 45. Again, the model f i t s many of the major features of the data. The largest discrepancy i s the reduced amplitude of the middle bump of the f i r s t two observed residuals. This may be due to the presence of another mode of o s c i l l a t i o n . 88 2% OF CONT i u I 4548 4552 4556 Wavelength i n angstroms Figure 44 Five model absorption lines versus observed lines from the A p r i l 19/20, 1982 data set. 89 {H5 % OF CONT 4 5 4 8 4552 4 5 5 6 Wavelength i n angstroms Figure 45 Model residuals plotted versus residuals from the A p r i l 19/20, 1982 data set. The phase difference between each model residual corresponds to the time difference between each observed residual. The model phase increments were calculated as (900 sees)' (2.54X10-1* rad s - 1 ) , producing an observed acceleration of the model absorption features of about 0.0055 km s - 1 (see Table VI). 90 Only a single mode, with 1=8 and m=-8, was used to produce the model l i n e s . I think i t unlikely that the addition of an 1=16 mode, which would produce a considerably more complicated p r o f i l e , would s i g n i f i c a n t l y improve the f i t . However, only a model that does incorporate this mode can provide the information necessary to decide this for certain. To summarize, the main conclusions that can be drawn from this work are: 1) Spica i s o s c i l l a t i n g i n a nonradial mode of order 1=8, m=-8. Additional modes of different orders are probably present as we l l . 2) Model results indicate that the ve l o c i t y amplitude of the o s c i l l a t i o n s i s small, of order 5 to 10 kms - 1 at most. 3) Accelerations of the observed absorption features imply an o s c i l l a t i o n angular frequency of between 3X10 - 5 and 9X10~5 r a d s - 1 . This frequency gives a period between 20 and 50 hours. Only a very high order G-mode would have a period of thi s length. 4) Matching model li n e s to the observed data results i n an estimated V of 190 kms - 1. This value compares with an eq expected synchronous rotation rate of about 100 kms -1 91 REFERENCES 1. Cox, John P., 1980, Theory of S t e l l a r Pulsation. (Princeton: Princeton University Press). 2. Cox, John P., 1984, P.A.S.P., 96_, 577. 3. Dukes, R.D., 1974, AP. J., 192, 81. 3. Fahlman, G.C. and Glaspy, J.W., 1973, Astronomical Observations with Television-Type Sensors, J.W. Glaspy and G.A.H. Walker, eds. 4. Gray, David F., 1976, The Observation and Analysis of S t e l l a r Photospheres (New York: John Wiley and Sons, Inc.). 5. Herbison-Evans, D., Hanbury-Brown, R., Davis, J . and Al l e n , J.R., 1971, M.N.R.A.S., 151, 161. 6. Hutchings, J.B. and H i l l , G., 1977, Ap. J . , 213, 111. 7. Lesh, J.R. and Aizenman, M.L., 1978, Ann. Rev. Astron. Astrophys., 16_, 215. 8. Lomb, N.R., 1978, M.N.R.A.S., 189 , 325. 9. Mihahlas, 1979, M.N.R.A.S., 189, 671. 10. Moyles, K., 1982, M.Sc. Thesis, University of B r i t i s h Columbia. 11. Osaki, Y., 1975, Pub. Astr. Soc. Jpn., 27_, 237. 12. Petr i e , R.M., 1962, "The Determin. of Orb. Elements of Spica Binary" Astronomical Techniques, Ed. W.A. H i l t n e r , Stars and St e l l a r Systems, Vol I I . 13. Saio, H., 1981, Ap. J . , 244, 299. 14. Shobbrook, R.R. , Lomb, N.R. and Herbison-Evans, D., 1972 , M.N.R.A.S., 156, 165. 15. Shobbrook, R.R., Herbison-Evans, D., Johnston, I.D. and Lomb, N.R., 1969, M.N.R.A.S., 145, 131. 16. Smith, M.A., 1979, Nonradial and Nollinear S t e l l a r Pulsations, ed. H.H. H i l l ands E.Z. Dziembouski (New York: Sprlnger-Verlag). 17. Smith, M.A., 1985, Ap. J., In Press. 18. Stamford, P.A. and Watson, R.D., 1976, Proc. Astron. Soc. Austr., _3, 75. 92 19. Stamford, P.A. and Watson, R.D., 1980, Proc. Astron. Soc. Austr., 4_, 80. 20. Unno, W., Osaki, Y. and Shibashi, H., 1979, Nonradial O s c i l l a t i o n of Stars. (Tokyo: Univ. Of Tokyo Press). 21. Vogt, S.S. and Penrod, G.D., 1983, Ap. J . , 275, 661. 22. Walker, G.A.H., Yang, S. and Fahlman, G.G., 1979, Ap. J., 233, 199. 23. Walker, G.A.H., Moyles, K., Yang, S. and Fahlman, G.G., 1982, P.A.S.P., 94, 143. 24. Walker, G.A.H., Johnson, R., and Yang, S., 1983, 8 t h Symp. on Photoelectronic Image Devices Proc. i n "Advances In Electronics and Electron Physics", B.L. Morgan, ed. (London: Academic Press). 25. Yang, S., 1980, M.Sc. Thesis, University of B r i t i s h Columbia. 93 APPENDIX I L i s t i n g of the program used to convert o r b i t a l elements into r a d i a l v e l o c i t y . Program also calculates the corresponding p i x e l s h i f t s between secondary and primary spectral lines 94 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC 1d = FRAS Page 1 1 rccccc 2 C This program calculates a nonradial o s c i l l a t i o n 3 C l l n e p r o f l l e using the supplied mode, th i s model 4 C does i t s c a l c u l a t i o n s around the complete sphere,le 5 C 1t ignores the limb. Any surface area that i s 6 C not visible,however, does not contribute to the l i n e p r o f i l e 7 C At the present time (March 07/85) only 1=6,7,8 1s supported. 8 C See subroutine MODE for the av a i l a b l e values of m for each 1. 9 C The t i d a l e f f e c t of a secondary is included as an 1=2, m=-2 mode. 10 C 1 1 C 12 C 13 C Program MAIN reads the Input data, and feeds t h i s to the model 14 C 1n turn, the model Is run by Subroutine CONTROL 15 C 16 REAL*8 AMPL(3),PH10(3),SIG(3),LAMO,VRMAX,XRMAX.K,PMAX 17 REAL*8 SPECT0(500),SPECTF(500),LAM(500).DLAM,ML).OME,OMAX,RMAX 18 REAL*8 OBL.RPOL.BETA.PI,INC,BPHIO,GRAVO,SPHIO 19 REAL*8 TI0E(3),TPHS(3),SIGTD,STPHS,TSTEP.SIGAPP 20 INTEGER NSTEP.NDEL,NALF,ENDS(2).NEL,LINDEX(7),LMODE,MMODE,STEP 21 22 CHARACTER FNAME1*22 23 CHARACTER FNAME2*11 24 CHARACTER FNAME3* 12 25 CHARACTER FNAME4*12 26 CHARACTER FNAME5*9 27 CHARACTER IMGNUM*2 28 29 COMMON /E1/AMPL,PHI0/E2/SIG 30 COMMON /E3/ENDS,LAMO,NSTEP 31 COMMON /E4/MU,OME,OBL.RPOL,BETA,PI.INC,BPHIO,GRAVO 32 COMMON /E5/SPECT0.SPECTF,LAM 33 COMMON /E6/LM0DE,MMODE 34 COMMON /E7/DLAM 35 COMMON /TIDE/SIGTD,TIDE,TPHS 36 37 READ (15.*) VRMAX,VPHS.K.LMODE,MMODE,VTIDE,SIGTD 38 READ (15,*) (PHIO(m),M=1,3) 39 READ (15,*) (TPHS(M),M=1,3) 40 READ (15.*) LAMO,DLAM 41 READ (15.*) NSTEP 42 43 C Calculate wavelengths of the Input spectrum 44 NEL=NSTEP+1 45 DO 2 1=1,NEL 46 2 LAM(I ) = LAMO+(1-1)*DLAM 47 48 READ (15,•)(SPECTO(I),1=1,NEL) 49 READ .(15,*)MU,OME.OBL,RPOL.BETA,INC.BPHIO,NDEL,NALF 50 READ (15,*) STEP.TSTEP 51 52 c ************************ 53 C 54 C Vrmax = maximum radial v e l o c i t y of the wave (kms) 55 C Vphs = apparent phase v e l o c i t y of the mode In radians/sec 56 C k = f r a c t i o n of vrmax that the tangential components can assume 57 C Lmode = L value of the o s c i l l a t i o n (1=6,7.8 presently available) 58 C Mmode = M mode of the o s c i l l a t i o n 95 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 2 59 C Vtlde = Maximum ra d i a l v e l o c i t y of the t i d a l o s c i l l a t i o n 60 C (assumed 1=2,m=2) 61 C Slgtd = angular frequency of the t i d a l o s c i l l a t i o n 62 C PhIO = s t a r t i n g phase for the o s c i l l a t i o n 63 C Tphs = s t a r t i n g phase for the t i d a l o s c i l l a t i o n 64 C LamO » s t a r t i n g wavelength of the Input spectral l i n e 65 C Dlam = wavelength step of the Input spectral l i n e ( In nm) 66 C Mu = Gr a v i t a t i o n a l constant (G*M) 67 C Ome = Rotational v e l o c i t y in radians/sec 68 C Obi = oblateness 69 C RPOL = polar radius in meters 70 C BETA = 1imb darkening parameter 71 C Inc = i n c l i n a t i o n to the l i n e of sight 72 C BPHIO = azimuth of the observers' l i n e of sight 73 C NDEL = number of d e c l i n a t i o n steps over 90 degrees 74 C NALF = number of azimuthal steps over 360 degrees 75 c ************************** 76 77 78 79 C Calculate the i n t r i n s i c angular frequency s1g(1) of the mode 80 C This 1s Just the d i f f e r e n c e between the angular frequency Ome due 81 C to r o t a t i o n , and the apparent phase v e l o c i t y (Vphs), times the 82 C m mode. The frequency that I'm c a l l i n g the i n t r i n s i c frequency 83 C here is not the real fundamental frequency being excited as I havn't 84 C included the c o r i o l i s c o rrection. 85 86 C Apparent angular frequency. This dictates 87 C how quickly the l i n e p r o f i l e s change 88 89 Slgapp = vphs*mmode 90 91 S1g(1)=Sigapp-mmode*Ome ! I n t r i n s i c angular frequency 92 93 CALL MAXCAL(LMODE,MMODE,PMAX,QMAX,RMAX) 94 CALL PARAM(VRMAX,XRMAX,K.PMAX,QMAX,RMAX,AMPL.SIG,SIGTD,VTIDE,TIDE) 95 96 C Phase increases with time due to the apparent 97 C angular frequency of the o s c i l l a t i o n 98 99 SPHIO=TSTEP*Sigapp 10O 101 STPHS=TSTEP*SIGTD ! phase step for the t i d a l contribution 102 103 WRITE (6,*) STEP.SPHIO 104 105 106 107 C Open the required output f i l e s , a new f i l e for each i t e r a t i o n 108 109 110 DO 10 J=1.STEP 1 1 1 112 IMGNUM=CHAR(J/10+48)// 113 & CHAR(d-10*(d/10)+48) 1 14 115 FNAME1='[FRASER.MODELOUT]OUT'//IMGNUM//'.DAT' 116 0PEN(16,NAME = FNAME1,STATUS='NEW' ) 96 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CCid=FRAS Page 3 117 FNAME2='TRACE'//IMGNUM//'.DAT' 1 18 OPEN(17,NAME = FNAME2,STATUS*'NEW) 1 19 FNAME3='SMRTRA'//IMGNUM//' .DAT' 120 OPEN(18,NAME=FNAMES,STATUS='NEW) 121 FNAME4='BINTRA'//IMGNUM//'.DAT' 122 OPEN(19,NAME=FNAME4,STATUS1'NEW) 123 FNAME5='MAP'//IMGNUM//'.DAT' 124 OPEN(20,NAME=FNAME5,STATUS1'NEW) 125 126 DO 20 1 = 1,500 • 127 SPECTF(I)=0.ODOO 128 20 CONTINUE 129 130 WRITE (16.499) LMODE,MMODE 131 WRITE (16,498) XRMAX 132 • WRITE (16,503) 133 WRITE (16,501) (AMPL(M),M=1,3) 134 WRITE (16.504) 135 WRITE (16,501) (PHIO(M).M=1.3) 136 WRITE (16,505) 137 WRITE (16,501) (SIG(M).M-1.3) 138 WRITE (16,506)LAM0.DLAM,NSTEP 139 ENDS(1)=0 140 ENDS(2)=NSTEP 14 1 142 WRITE (16,509) 143 WRITE (16,510)(LAM( I ) ,SPECT0(I),1 = 1,NEL) 144 GRAVO=-MU/(RPOL*RPOL) 145 WRITE (16,511)MU,OME.OBL.RPOL,BETA,INC,BPHIO.GRAVO, 146 0 NDEL,NALF 147 PI=DAC0S(-1.OD+OO) 148 CALL CONTROL(NDEL.NALF.NSTEP+1) 149 150 CL0SE(16) 151 CL0SE(17) 152 CL0SEO8) 153 CL0SEO9) 154 CL0SE(2O) 155 156 C Increment the mode phase and the t i d a l phase. 157 C The phase Increment 158 C includes the increment due to the wave propagation, and the 159 C r o t a t i o n of the star (OME) 160 161 PHIO(1)=PHIO(1)+SPHIO 162 PHIO(2)=PHIO(2)+SPHIO 163 PHIO(3)=PHIO(3)+SPHIO 164 165 TPHS(1)=TPHS(1)+STPHS 166 TPHS(2)=TPHS(2)+STPHS 167 TPHS(3)=TPHS(3)+STPHS 168 169 C 170 171 172 10 CONTINUE 173 174 498 FORMAT(' MAXIMUM EQUATORIAL EXCURSION : ',012.6) 97 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 4 175 499 FORMAT(' MODE USED FOR THIS MODEL: L='.I3,' M='.I3) 176 501 FORMATC RADIAL WAVES=',D12.6// 177 0' DECLINATION WAVES=',D12.6// 178 0' RIGHT ASCENSION WAVE='D12.6//) 179 503 FORMAT('1 AMPLITUDES OF THE WAVE COMPONENTS 180 0(METERS): './/) 181 504 FORMAT('1 PHASES OF THE WAVE COMPONENTS: ',//) 182 505 FORMATC1 ANGULAR FREQUENCY OF THE WAVE COMPONENTS 183 Q(RAD/S) './/) 184 506 FORMAT('1 SPECTRUM:'//' STARTING WAVELENGTH: 185 Q.F10.4,' NM'/' STEP WAVELENGTH: '.F10.4,' NM'/ 186 Q' NUMBER OF STEPS: ',15) 187 509 FORMAT(/' WAVELENGTH-INTENSITY PAIRS:'/) 188 510 F0RMAT(4(F13.6,'='.F13.6)) 189 511 FORMAT('1 MISCELLANEOUS:'//' GRAVITATIONAL PARA' 190 Q, 'METER OF THE STAR: ',4X,D 13.7, ' M3/S2'/' ROTATION' 191 Q.'AL VELOCITY :' ,20X.D13.7, ' RAD/S'/' OBLATENESS: ' 192 Q.29X.D13.7/' POLAR RADIUS: ' , 27X,D 13.7,' M'/ 193 Q' LIMB DARKENING PARAMETER: ' . 11X,F10.4/ 194 Q' INCLINATION TO LINE OF SIGHT:',7X,F10.4,' RAD'/ 195 Q' GLOBAL PHASE:',23X,F10.4/' POLAR GRAVITY:',22X 196 Q.F10.4,' M/S2'//' NUMBER OF DECLINATION STEPS:',8X, 197 QI5,' OVER 90 DEGREES'/' NUMBER OF RIGHT ASCENSION' 198 Q,' STEPS:',4X.15,' OVER 360 DEGREES') 199 STOP 200 END 201 202 C 203 C Subroutine CONTROL runs the other subroutines 204 C 205 SUBROUTINE CONTROL(NDEL,NALF,NEL) 206 C 207 INTEGER NDEL,NALF,NEL,ENDS(2) 208 INTEGER NALD(90),CPTR.SDE,SDAL.OLCPTR 209 REAL *8 MATCO(3,4),MATPO(3,4),MATC1(3,4),XI(3 . 4) 210 REAL*8 VOBPO) , V0BC(3) ,RAD(2) . VS(3) , VC(3) . VP(3) , VXI (3) 211 REAL*8 Z,G,S,DE,DDE,BIDON(5),PI,INC,BPHIO,VXIP(3) ,VEL 212 REAL*8 GRAVO,DAL(90),R,KS 213 COMMON /E3/ENDS,LAMO,NSTEP 214 COMMON /E4/BID0N,PI,INC,BPHIO,GRAVO 215 COMMON /E7/DLAM 2 16 217 C 218 C Calculate the number of RA steps per declination interval using 219 C the r e l a t i o n #of steps=# of RA steps Input*cos((1-1)*dec) 220 C and then c a l c u l a t e the number of radians per right ascension steps 221 C 222 CPTR=0 223 SDE=1 224 DDE = PI/(2.*FL0AT(NDEL)) 225 DO 1 1=1.NDEL 226 NALD(I ) = INT(FLOAT(NALF)*COS(FLOAT(I - 1)*SNGL(DDE))) 227 1 DAL(I)=2.D+0O*PI/DBLE(FL0AT(NALD(I))) 228 WRITE (16,100)(NALD(I),DAL(I),I,1=1,NDEL) 229 WRITE (16.101) 230 VOBP(1)=1.D+OO 231 V0BP(2)=INC 232 V0BP(3)=BPHIO 98 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 5 233 C 234 C VOBP contains the polar coordinates of the o r i g i n 235 C VOBC contains the c a r t e s i a n coordinate of the o r i g i n 236 C 237 CALL POLCA(VOBP.VOBC, 1) 238 DE=PI/2.0+00 239 C 240 1007 DO 1000 IDE= 1,NDEL 241 0LCPTR=CPTR 242 C 243 C Ca l c u l a t e s t e l l a r radius at the top and bottom of the surface 244 C being considered 245 C 246 SDAL=1 247 CALL RADIUS(DE.RAD(1+(1-SDE)/2)) 248 CALL RADIUS(0E-DDE,RAD(2-(1-SDE)/2)) 249 C 250 C Put the polar coordinates of the undisturbed four corners of the 251 C uni t surface under consideration into MatpO 252 c 253 DO 2 1=1,2 254 MATP0(1.I)=RAD(1) 255 MATPO(1,I+2)=RAD(2) 256 MATP0(2,I+1-SDE)=DE 257 MATP0(2.I+1+SDE)=DE-DDE 258 MATP0(3.3*I-2)=BPHI0 259 2 MATPOO. I-M)=BPHIO+DAL(IDE) 260 c 261 c Calculate the I n i t i a l XI values for the f i r s t unit area at t h i s 262 c d e c l i n a t i o n XI contains the displacement of the four corners 263 c 264 1003 CALL DECIDE(MATPO.XI.O) 265 CALL POLCA(MATPO,MATCO,4) 266 C 267 C Add the disturbance due to NR wave to i n i t i a l undisturbed coord. 268 C 269 DO 3 1=1,4 270 DO 3 d=1.3 271 3 MATC1(J,I)=MATC0(d,I)+XI(U.I) 272 IIAL=NALD(IDE)/2 273 KS=DAL(IDE)*DTAN(DDE/2.D+0O)/(D0E*DSIN(DAL(IDE)/2.D+OO)) 274 C 275 C Do 10O1 from 1 to the H of ri g h t ascension steps for t h i s dec 276 C 277 DO 1001 IAL=1.IIAL 278 CALL SURF(MATC1,S,VS,VC) 279 S=S*KS 280 CALL DOTC(VS.VOBC.R) 281 C 282 C R= the p r o j e c t i o n of the average surface perpendicular onto the 1 283 C sight to the observer 284 C S= the surface area of the unit area 285 C KS= the c o r r e c t i o n f or departures from a spherical surface 286 C 287 C 288 IF (SNGL(R).LE.0.) GO TO 10O1 289 C 1e If we can't see the surface (R less than 0) go to 1001 290 C 99 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CCid=FRAS Page S 291 CALL CARPO(VC,VP,1) 292 CALL SPEED(VP.VXI,VXIP) 293 CALL GRAV(VP.G.VXIP) 294 CALL GDARK(G.S) 295 CALL LDARK(R.S) 296 CALL DOTC(VXI,VOBC.Z) 297 Z = -Z 298 VEL=Z 299 CALL HISTO(VEL.O) 300 CALL SPEC(Z,(R*S)) 301 CPTR=CPTR+1 302 WRITE(16. 110)CPTR,IDE.IAL,ENDS,S.Z,G,R,VP(2).VP(3) 303 CALL MAPER(VP(2),VP(3).X.Y,INC) 304 WRITE(20.666) X.Y.Z 305 306 DO 8 1=1,3 307 MATC1(I.2-SDAL)=MATC1(I.3-SDAL) 308 MATC1(1.3+SDAL)=MATC1(1.2+SDAL) 309 MATCO(I.2-SDAL)=MATC0(I.3-S0AL) 310 8 MATCO(I.3+SDAL)=MATCO(I,2+SOAL) 311 DO 7 1=1,4 312 7 MATPO(3,I)=MATPO(3.I)+DAL(IDE)*SDAL 313 CALL DECIDE(MATPO,XI,SDAL ) 314 DO 4 1=1,2 315 >J=(2-SDAL)*I + (3*SDAL-1 )/2 316 MATCO(1,d)=MATPO(1,J)*DSIN(MATPO(2.J))*DCOS(MATPO(3 . J)) 317 MATC0(2.J)=MATPO(1,J)*DSIN(MATPO(2,d) )*DSIN(MATPO(3,J)) 318 DO 4 K=1,3 319 4 MATC1(K,J)=MATCO(K,j)+XI(K,J) 320 10O1 CONTINUE 321 C 322 C Having c a l c u l a t e d the co n t r i b u t i o n from a declination s t r i p , we e x i t the 323 C the 1001 loop. If the d e c l i n a t i o n Just calculated made no co n t r i b u t i o n 324 C to the spectrum, OLCPTR w i l l equal CPTR. In this case we assume we have 325 C reached a northern, or southern limb. (If the southern limb, then SDE 326 C w i l l equal -1 and we e x i t ) 327 C 328 IF(OLCPTR.EO.CPTR) GO TO 1004 329 IF (SDAL.EQ.-1) GO TO 1005 330 SDAL=-1 331 WRITE (16,105) ! Completed a half d e c l i n a t i o n s t r i p 332 DO 6 1=1,2 333 MATPO(3.3*I-2)=BPHIO-DAL(IDE) 334 6 MATP0(3,1+1)=BPHIO 335 GO TO 1003 336 1005 DE=PI/2.D+OO-DBLE(FLOAT(IDE))*DDE 337 WRITE (16.106) ! Going to the next d e c l i n a t i o n 338 1000 CONTINUE 339 1004 IF (SDE.EQ.-1) GO TO 1006 340 WRITE (16.107) ! Completed the northern hemisphere 341 SDE=-1 342 DDE=-DDE 343 DE=PI/2.D+00 344 GO TO 1007 345 1006 CONTINUE 346 WRITE (16.108) ! Completed the whole star 347 WRITE (16.109) CPTR 348 CALL TRACE(ENDS(2)+1-ENDS(1)) 100 L i s t i n g Of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CCi d-FRAS Page 7 349 VEL=Z 350 CALL HIST0(VEL,1) 351 352 100 FORMATC1 NUMBER OF RIGHT ASCENSION STEPS AND ' 353 0.'LENGTH OF THE STEPS VS DECLINATION INDEX'/(I10 354 Q.F15.10,110)) 355 101 FORMATC1 CPTR IDE IAL ENDS',11X,'S'.16X,'Z',15X 356 0, ' G ' , 15X , ' R ' 15X , 'DEC , 15X, ' RA ' / ) 357 104 FORMATC COMPLETED A TURN WITHOUT CROSSING LIMB') 358 105 FORMAT( ' COMPLETED A HALF DECLINATION STRIP') 359 106 FORMATC . GOING TO THE NEXT DECLINATION') 360 107 FORMAT(/' COMPLETED NORTHERN HEMISPHERE'/) 361 108 FORMAT(/' COMPLETED WHOLE STAR') 362 109 FORMAT(///' CUMULATED UNIT SURFACES: '.16) 363 110 F0RMAT(5I6,6D16.6) 364 666 F0RMAT(2F6.3,F15.2) 365 RETURN 366 END 367 368 subroutine apconv (vector 1,npt1,vector2,npt2,vector3,nout) 369 370 C This routine convolves two vectors using the array processor 371 C Courtesy of G. Grieve 372 C C a l l i n g sequence: 373 374 real vector 1(npt1) ! In :operand vector 375 Integer npt1 !in :U of elements in vectorl 376 real vector2(npt2) !1n operator vector 377 integer npt2 !In :# elements in vector2 378 real vector3(npt 1 -npt2+1) ! out -.reultant vector 379 integer nout ! out :H of elements out 380 381 integer v1_add !address of vectorl 382 integer v2_add !address of vector2 383 integer v3_add ! address of vec.tor3 384 385 c a l l a p i n i t (0.0.1) 386 387 388 v1_add = O 389 v2_add = v1_add + npt1 390 v3_add = v2_add + npt2 391 nout = npt1 - npt2 + 1 392 c a l l apput (vector 1.v1_add.npt1,2) 393 c a l l apput (vector2,v2_add.npt2,2) 394 cal1 apwd 395 396 c a l l conv (v1_add,1, (start address of v1 397 & v3_add-1.-1, !end address of v2 398 & v3_add,1,nout,npt2) 399 400 cal1 apwr 401 402 c a l l apget (vector3,v3_add,nout,2) 403 c a l l apwd 404 405 c a l l aprlse 406 / 101 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 8 407 408 return 409 end 410 , 411 412 SUBROUTINE ATIDE(DE.AL,DXI,PH) 413 C 414 C Calculate the displacements at the corners of the input unit 415 C area for the 1=2, m = + 2 mode, being used to simulate the t i d a l d i s t o r t i o n 416 C 417 ' 418 REAL*8 DE,AL,DXI(3),PH,PI,BID0N(5).DCONST,M1 419 REAL*B TIDE(3),TPHS(3),X,Y,SIGTD 420 COMMON /TIDE/SIGTD,TIDE.TPHS 421 C L=2,M=+2 422 423 X=DCOS(DE) 424 Y=DSIN(DE) 425 426 DXI( 1 ) = -TIDE(1)*DSIN(TPHS(1) + 2*AL-PH*PI)* 427 & 3.DOO*Y**2 428 DXI(2)=-TIDE(2)*DSIN(TPHS(2)+2*AL-PH*PI)* 429 & 6.D00*X*Y 430 431 DXI(3)=-TIDE(3)*2*DC0S(TPHS(3)+M1*AL-PH*PI)* 432 & (-6.D00*Y) 433 RETURN 434 END 435 436 437 C 438 SUBROUTINE BINDAT(DATA.BIN,N) 439 C 440 C This routine averages f i v e p i x e l s together, and outputs the result 441 C as a rebinned trace. This accounts for the fact that the model 442 C c a l c u l a t e s the l i n e at 0.03 angstrom spacings while the data is 443 C taken at 0.15 angstrom spacings. Note that the model spacing is an 444 C input parameter and can be changed. 445 446 REAL*4 DATA(512),BIN(10O) 447 INTEGER N.N1 448 N1=60 449 DO 10 1=1.100 450 - BIN(i)=1.00 451 10 continue 452 453 11=1 454 DO 20 1=10,70 455 BIN(I) = ((DATA(I 1)+DATA(I 1 + 1)+DATA(I 1+2)+DATA(I 1 + 3) + 456 S DATA(I 1 + 4))/5.0) 457 11=11+5 458 20 CONTINUE 459 RETURN 460 END 461 462 C 463 SUBROUTINE CARPO(MC.MP,J) 464 C 102 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 9 465 C Convert from cartesian to polar coordinates 466 C 467 REAL*8 MC(3.d),MP(3,d) 468 DO 1 K=1.d 469 MP(1,K)=0.0+00 470 DO 2 1=1.3 471 2 MP( 1 ,K)=MP(1,K)+MC(I,K)*MC(I,K) 472 MP( 1 ,K)=DSQRT(MP(1.K)) 473 MP(2.K)=DAC0S(MC(3,K)/MP(1,K)) 474 1 MP(3,K)=DATAN2(MC(2,K),MC(1,K)) 475 RETURN 476 END 477 478 c 479 SUBROUTINE CROSSC(VA.VB.VC ) 480 c 481 c Perform a cross product of the Input vectors VA and VB 482 c VC 1s returned with the r e s u l t 483 c 484 REAL*8 VA(3),VB(3),VC(3) 485 INTEGER IN(4)/2,3,1,2/ 486 DO 1 1=1.3 487 1 VC(I)=VA(IN(I))*VB(IN(I+1))-VA(IN(I+1))*VB(IN(I)) 488 RETURN 489 END 490 491 492 c 493 SUBROUTINE DECIDE(MATPO.XI.SDAL ) 494 c 495 c 496 c Decide i f th i s is the f i r s t unit area in a declinat 497 c (SDAL=0),or a subsequent area to the west (SDAL=1) 498 c (SDAL=-1) 499 c 500 INTEGER SDAL,I 1.12.13,IL 501 REAL*8 MATP0(3,4).XI(3,4) ,MEM(2) ,DE, AL.DXIO) 502 INTEGER LINDEX(9) 503 11 = 1 504 12 = 4 505 13= 1 506 507 IF (SDAL) 4,8,6 508 8 DO 9 I=1,4 509 DO 9 d=1 ,3 510 9 XI(J,I)=0.D+OO 51 1 GO TO 5 512 6 DO 1 1=2.3 513 DO 1 J=1.3 514 XI(d.3*I-5)=XI(d,I) 515 1 XI(d,I)=0.D+OO 516 I 1=2 517 12 = 3 518 GO TO 5 519 4 DO 7 1=2,3 520 DO 7 d=1,3 521 XI(d,I )=XI(d,3*1-5) 522 7 XI(d.3*I-5)=O.D+00 103 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 10 523 13 = 3 524 5 DO 2 1=11,12,13 525 DE=MATP0(2.1) 526 AL=MATPO(3 , I) 527 C 528 C C a l l MODE to c a l c u l a t e the displacements at the corners 529 C unit area 530 C 531 CALL M0DE(DE,AL,DXI.0.D+0O) 532 DO 3 JM.3 533 XI(J,I)=XI(J,I)+DXI(J) 534 3 CONTINUE 535 536 CALL ATIDE(DE,AL.DXI.O.OD+OO) 537 , DO 33 J*1.3 538 XI(J,I)=XI(d,I)+DXI(J) 539 33 CONTINUE 540 C 54 1 C Convert the displacements from an R.Theta, Ph1 system to 542 C the usual X,Y,Z,system 543 C 544 MEM(1)=DCOS(AL)*(XI(1,I )*DSIN(DE)+XI(2.I)*DCOS(DE) ) 545 0-DSIN(AL)*XI(3,I) 546 MEM(2)=DSIN(AL)*(XI(1,I ) *DSIN(DE)+XI(2,I)*DCOS(DE)) 547 Q+DCOS(AL)*XI(3,I) 548 XI(3,I)=DCOS(DE)*XI(1,I)-DSIN(DE)*XI(2.I) 549 XI(2,I)=MEM(2) 550 2 XI(1,I)=MEM(1) 551 RETURN 552 END 553 554 C 555 SUBROUTINE DOTC(VA,VB,R) 556 C 557 C Do a dot product on VA and VB, R contains the re s u l t 558 C 559 REAL*8 VA(3),VB(3).R 560 R=0.D+0O 561 DO 1 1=1,3 562 1 R=R+VA(I)*VB(I) 563 RETURN 564 END 565 566 C 567 SUBROUTINE GDARK(G.S) 568 C 569 c Do g r a v i t y darkening 570 c 571 REAL*8 G,S,GRAVO.BID0N(8) 572 COMMON /E4/BID0N.GRAVO 573 S=DABS(S*G/GRAVO) 574 RETURN, 575 END 576 577 c 578 SUBROUTINE GRAV(VP,G,VXIP) 579 c 580 c C a l c u l a t e local g r a v i t y 104 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 11 581 C 582 INTEGER IL 583 RE AL *8 SIG( 3 ), VP( 3 ) , G, Mil, VXIP ( 3 ) . DX I ( 3 ) 584 REAL*8 SIGTD,DX11(3) 585 COMMON /E2/SIG/E4/MU 586 COMMON /TIDE/SIGTD,TIDE,TPHS 587 G=O.D+00 588 589 CALL M0DE(VP(2),VP(3),DXI.O.D+OO) 590 CALL ATIDE(VP(2),VP(3),DXI 1 .O.D+OO) 591 592 , G=G-(DXI( 1 )*SIG(1)*SIG(1)+DXI1 ( 1)*SIGTD*SIGTD) 593 G=G-MU/(VP(1)*VP(1))+VXIP( 2)*VXIP(2)/VP( 1 ) + VXIP(3)* 594 QVXIP(3)/(VP(1)*DSIN(VP(2))) 595 RETURN 596 END 597 598 SUBROUTINE HISTO(VEL.FLAG) 599 C This subroutine calculates a histogram of number of unit surfaces versus 600 C v e l o c i t y . 601 INTEGER BIN(40),NBIN(40).FLAG 602 REAL*8 VEL 603 REALM Z1 604 VEL=VEL/10O0.D0O 605 Z1=SNGL(VEL) 606 IF(FLAG.EO.1)G0 TO 90 607 IF(Z1.LT.O.OEOO) GO TO 20 608 I=IFIX(1.0+Z1/5.0) 609 BIN(I)=BIN(I )+ 1 610 RETURN 611 20 Z1=ABS(Z1) 612 I=IFIX(1.0+Z1/5.0) 613 NBIN(I)=NBIN( I)+1 614 RETURN 615 90 N=40 616 WRITE(6,500) N 617 WRITE(6.501) (BIN(I),I=1.40) 618 WRITE(6,500) N 619 WRITE(6.501) (NBIN(I).I=1,40) 620 DO 30 1=1.40 621 BIN(I)=0 622 NBIN(I)=0 623 30 CONTINUE 624 625 500 FORMAT(13) 626 501 F0RMAT(5I5) 627 RETURN 628 END 629 630 631 632 633 634 635 636 637 638 SUBROUTINE LDARK(R.S) 105 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 12 639 C 640 C Do 1Imb darkening 64 1 C 642 REAL*8 R.S.BETA,BID0N(4) 643 COMMON /E4/BID0N.BETA 644 S = S*(1.D+00-BETA*(1 D+OO-R)) 645 RETURN 646 END 647 648 SUBROUTINE MAPER(NPD.LONG,X.Y.INC) 649 C Subroutine maper converts input l a t i t u d e (declination) and longitude 650 C to X and Y coordinates. Please note that latitude in the main program 651 C is measured from the north pole rather than from the equator, as is 652 C usual. This is usually c a l l e d the North Polar Distance (1e. NPD) 653 C Thus we make 1 atitude=(90-NPD) . S i m i l a r l y , the longitude is measured 654 C from the observers meridian, we need longltude=(90+LONG). The input 655 C i n c l i n a t i o n of the star is the angle from the stars pole to the 656 C the observer defined equator. We need the displacement of the s t a r ' s 657 C pole from the observer defined pole. This displacement is given by 658 c E with E=(90-INC) 659 660 661 REAL*8 LAT,LONG,XD,YD,NPD,PI,E,INC 662 REAL*4 X.Y 663 664 PI=3.1415927D00 665 PI2=PI/2.DOO 666 667 L0NG=(L0NG-PI2) 668 LAT=(PI2-NPD) 669 E=(PI2-INC) 670 671 X1=DC0S(LAT)*DC0S(L0NG) 672 Y1=DC0S(LAT)*DSIN(E)*DSIN(LONG)+DC0S(E)*DSIN(LAT) 673 674 X=(1.0+SNGL(X1))/2.0 675 Y = (1.0+SNGL(Y1))/2 .0 676 RETURN 677 END 678 679 SUBROUTINE MAXCAL(L.M,PMAX,OMAX,RMAX) 680 c 681 c Calculates the maximum value of the Legendre polynomial for 682 c for various modes. At the moment only 1=8, m=8,7,6,5,4,3,2 or 683 c 1=6, m=6,5,4,3.2 .or 1=7, m=7.6,5.4,3 , are available. 684 c 685 INTEGER L,M,N 686 REAL*8 P(400),0(400),R(400),X,Y 687 REAL*8 PMAX,OMAX,RMAX,STEP,PDEG,RDEG.ODEG.DE,PI 688 PI=DAC0S(-1.000) 689 STEP=PI/200.DOO 690 PMAX=O.ODOO 691 QMAX=O.ODOO 692 RMAX=0.0D00 693 DE=0.0000 694 M1 =DFL0AT(M) 695 N=IABS(M) 696 DO 100 1=1,200 106 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC i d = FRAS Page 13 G97 IF((L+1 )-8)6.7,8 698 C L = 8, M = + OR - 8 699 8 IF(N+1-8)86,87.88 70O 701 88 P(I)=2027025.DOO/128.DOO*(35.DOO+DC0S(8.ODOO*DE)-8.DOO* 702 & DC0S(6.DOO*DE)+28.DOO*DC0S(4.DOO*DE)-56.DOO* 703 & DC0S(2.D00*DE)) 704 705 0(I)=2027025.DOO/128.DOO*(-8.DOO*DSIN(8.DOO*DE) + 48 .DOO* 706 & DSIN(6.D0O*DE)-112.DOO*DSIN(4.DOO*DE)+112.DOO* 707 & D5IN(2.D0O*DE)) 708 709 R(I)=M1*2027025.DOO*(DSIN(DE)**7+6.DOO'DSIN(OE)) 710 GO TO 90 71 1 C L = 8, M=+ OR -7 712 713 87 X=DCOS(DE) 714 Y=DSIN(DE) 715 P(I)=2027025.DOO*Y*(-X**7+3 DOO*X**5-3.DOO*X**3+X) 716 717 0(I)=2027025.D00*((Y**2)*(7.D00*X**6-15.D00*X**4+9.DOO*X**2-718 & 1.DOO)+(-X**8+3.DOO*X**6-3.DOO*X**4+X**2)) 719 720 R(I)=M1*2027025.DOO*(-X**7+3.DOO*X**5-3.000*X**3+X) 721 722 GO TO 90 723 86 IF(N+1-6)84.85,886 724 C L=8 , M=+ OR -6 725 886 X=DCOS(DE) 726 Y=DSIN(DE) 727 P( I ) = 135135.DOO/2.DOO*(- 15.DOO*X**8+46.DOO*X**6-48 . D00*X**4 728 & +18.DOO*X**2-1.DOO) 729 730 0( I) = 135135.DOO/2 .D0O*Y*(120.DOO*X**7-276.DOO*X**5+192.DOO* 731 & X**3-36.DOO*X) 732 733 R(I)=M1 * 135135.DOO/2.DOO*(14.DOO*Y**5-15.DOO*Y**7) 734 735 GO TO 90 736 737 C L = 8, M=+ OR -5 738 85 X=DCOS(DE) 739 Y=DSINI(DE) 740 P(I)=135135.DOO/2.D0O*Y*(5.DOO*X**7-11.DOO*X**5+7 .DOO*X**3 741 & -x) 742 743 0( I )=135135.DOO/2.DOO*((5.DOO*X**8-11 .DOO*X**6+7.DOO*X**4-744 & X**2)-Y**2*(35.DOO*X**6-55 .DOO*X**4+21.DOO*X**2-1.DOO)) 745 746 R(I)=M1 * 135135.DOO/2.DOO*(5.D00*X**7-11.DOO*X**5+7.D0O*X**3-X) 747 748 GO TO 90 749 750 84 IF(N+1-4)82.83,884 751 C L = 8. M=+ OR -4 752 884 X=DCOS(DE) 753 Y=DSIN(DE) 754 107 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 14 755 P(I)=10395.D00/8.000*(65.D00*X**8-156 ,D00*X**6+118.DOO*X**4 756 & -28,000*X**2+1.000) 757 758 0(I)=-10395.D0O*Y*(65 D00*X**7-117.DOO*X**5 + 59.D00*X**3 759 & -7.D00*X) 760 761 R(I ) =M1*10395.DOO/8.DOO*(65.DOO*Y**7-104.DOO*Y**5+40.DOO*Y**3) 762 763 GO TO 90 764 83 X=DCOS(DE) 765 Y=DSIN(DE) 766 767 P(I)=3465.DOO/8.000*Y*(-39.DOO*X**7+65.DOO*X**5-768 & 29.DOO*X**3+3.DOO*X) 769 770 0(I)=3465.DOO/8.D00*((-39.DOO*X*+8 + 65.DOO*X**6-29 . DOO*X**4 + 771 & 3.DOO*X**2)+Y**2*(273.DOO*X**6-325.DOO*X**4+87.DOO*X**2 772 & -3.000)) 773 774 R(I)=M1*3465.DOO/8.DOO*Y*(-39.DOO*X**7+65.DOO*X**5-775 & 29.DOO*X**3+3.DOO*Xj 77ft 777 GO TO 90 778 779 C L=8. M=+ OR -2 780 82 X=DCOS(DE) 781 Y=D5IN(DE) 782 783 P(I)=315.DOO/16.DOO*(-143.DOO*X**8+286.DOO*X**6-176.DOO* 784 & X**4+34.DOO*X**2-1.DOO) 785 786 0(1 ) = 315.DOO/4.DOO*(-Y)*(286.DOO*X**7 + 429.DOO*X**5-176.DOO* 787 & X**3+17.DOO*X) 788 789 R(I)=M1*315.DOO/16.DOO*(- 143.DOO*Y**7+286.DOO*Y**5-176.DOO* 790 & Y**3+32.DOO*Y) 791 792 GO TO 90 793 C L=7, M=+ OR - 7 794 7 IF(N+1-7)75,76,77 795 796 77 X=DC0S(DE) 797 Y=DSIN(DE) 798 P(I ) = 1 35 135.DOO*Y*(-X**6+3.DOO*X* *4-3.D00*X* *2 +1.DOO) 799 800 0(1)=135135.DOO*(-(Y**2)*(-6.DOO*X**5+12.DOO*X**3-6.DOO*X)+ 801 & (-X**7+3.DOO*X**5-3.DOO*X**3+X)) 802 803 ' R(I)=M1*135135.D00*(-X**6+3.D00*X**4-3.DOO*X**2+1.DOO) 804 GO TO 90 805 806 C L=7, M=+ OR -6 807 808 76 X=DCOS(DE) 809 Y=DSIN(DE) 810 P(I)=135135.DOO*(-X**7+3.DOO*X* *5-3.DOO*X**3+X) 811 812 0(I)=135135.DOO*(-Y)*(-7.DOO*X**6+15.DOO*X**4-9.D00*X**2+ 108 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CCid=FRAS Page 15 813 & 1.DOO) 814 815 R(I)=M1*135135.DOO*Y*(X**5-2.DOO*X**3+X) 816 817 GO TO 90 818 75 IF(N+1-5)73.74.775 819 820 C L=7, M=+ OR - 5 821 822 775 X=DCOS(DE) 823 Y =DSIN(DE) 824 P( I)=10395.000/2.DOO*(Y)*(13.DOO*X**6-27.DO0*X**4+15.DOO* 825 & X**2-1.D00) 826 0( I) = 10395.DOO/2.DOO*(-(Y**2)*(78.D00*X**5-108.DOO*X**3 + 827 & 30.D00*X)+(13.DOO*X**7-27.DOO*X**5+15.D00*X**3-X)) 828 829 R(I)=M1*10395 DOO/2.DOO*(13.DOO*X**6-27.DOO*X**4+15.DOO* 830 & X**2-1.DOO) 831 B32 GO TO 90 833 834 C L=7, M=+ OR - 4 835 836 74 X=DCOS(DE) 837 Y=DSIN(DE) 838 P(I ) = 3465.DOO/2.DOO*(13.DOO*X**7-29.DOO*X**5+19.DOO* 839 & X**3-3.DOO*X) 840 841 0( I)=3465.DOO/2.DOO*(-Y)*(91.DOO*X**6-145.DOO*X**4 + 57.DOO* 842 & X**2-3.DOO) 843 844 R( I)=M1«3465.DOO/2.000*(Y)*(- 13.DOO*X**5+16.DOO*X**3-845 & 3.D00*X) 846 847 GO TO 90 848 849 73 IF (N+1-3)71.72,773 850 851 C L=7, M=+ OR -3 852 853 773 X=DCOS(DE) 854 Y=OSIN(OE) 855 P( I)=315.DOO/8.DOO*(Y)*(-143.DOO*X**6 + 209.DOO*X**4-69.DOO* 856 & X**2+3.DOO) 857 858 0( I )=315.DOO/8.DOO*(-(Y**2)*(-858.DOO*X* *5+836.DOO*X* *3-859 & 138.DOO*X)+(-143.0OO*X**7+209.DOO*X**5-69.DOO*X**3+ 860 & 3.DOO*X)) 861 862 R(I)=M1*315.DOO/8 DOO*(- 143.00O*X**6 + 2O9.DOO*X**4-863 & 69 0OO*X**2+3.000) 864 865 GO TO 90 866 867 C L=7, M=+ OR -2 868 869 72 X=DCOS(DE) 870 Y=DSIN(DE) 109 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CCid=FRAS Page 16 871 GO TO 90 872 C L = 7 . M=+ OR -1 873 874 71 X=DCOS(DE) 875 Y=DSIN(DE) 876 GO TO 90 877 878 879 6 IF(N+1-6)664.65,66 880 C L =6. M=+ OR -6 881 882 66 P(I) = 1 1.000*945.DOO/32.DOO*(-DCOS(6.DOO*DE)+6.DOO* 883 & DC0S(4.DOO*DE)-15.DO0*DC0S(2.DOO*DE)+1O.DOO) 884 885 0(I )=66.DOO*945.DOO/32.D00*(DSIN(6.DOO*DE)-4.DOO*DSIN(4. DOO* 886 & +5.DOO*DSIN(2.DOO*DE)) 887 888 R(I)=M1 * 11.000*945.DOO/16.DOO*(DSIN(5.DOO*DE)-5.DOO* 889 & DSIN(3.DOO*DE)+10.DOO*DSIN(DE)) 890 GO TO 90 891 C L =6, M=+ OR -5 892 893 894 895 65 P(I)=11.DOO*945.DOO/32.DOO*(DSIN(6.D0O*DE)+4.DOO*DSIN(4. DOO*l 896 & +5.DOO*DSIN(2.D0O*DE)) 897 898 0(I)=11.DOO*945.DOO/32.DOO*(6.DOO*DCOS(6.DOO*DE)+16.DOO 899 & *DCOS(4.DOO*DE)+10.DOO*DC0S(2.DOO*DE)) 900 901 R(I)=M1* 11.000*945.DOO/16.DOO*(DCOS(5.DOO*DE)-3.DOO* 902 & DCOS(3.DOO*DE)+2.DOO*DC0S(DE)) 903 GO TO 90 904 905 664 IF(N+1-4)62,63,64 906 C L = 6, M=+ OR -4 907 64 P(I)=9.D00*105.D00/64.D00*(11.DOO*DCOS(6.DOO*DE)-26.DOO 908 & *DCOS(4.DOO*DE)+5.DOO*DCOS(2.DOO'DE)+10.DOO) 909 910 0(I)=9.DOO*105.000/64.DOO*(-66.DOO*DSIN(6.DOO*DE)+104.DOO 91 1 & *DSIN(4.DOO*DE)-10.DOO*DSIN(2.DOO*DE)) 912 913 R(I)=M1*9.DOO*105.DOO/32.DOO*(- 11.D0O*DSIN(5.DOO*DE)+15. DOO* 914 & DSIN(3.DOO*DE)+1O.DOO*DSIN(0E)) 915 GO TO 90 916 C L = =6, M=+ OR -3 917 63 P(I) = 315.DOO/64.DOO*(- 11 .DOO*D5IN(6.DOO*DE)+12.DOO* 918 & DSIN(4.DOO*DE)+9.DOO*DSIN(2.DOO*DE)) 919 920 0(I) = 1890.DOO/64.DOO*(- 11.DO0*DC0S(6.DOO*DE) + 8.DOO* 92 1 & DC0S(4 DOO*DE)+3.DOO*DC0S(2 DOO'DE)) 922 923 R(I)=M1*315.DOO/32.DOO*(55.DOO*DC0S(0F)+DC0S(3.0OO*DE)-924 & 11.DOO*DC0S(5.D0O*DE)) 925 GO TO 90 926 C L = 6, M=+ OR -2 927 62 P(I)=105.DOO/256.DOO*(10.D00+17.DOO*DC0S(2.DOO*DE)+6.DOO* 928 & DCOS(4.DOO*DE)-33.DOO*DC0S(6.DOO*DE)) 110 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 17 929 930 0(1)=105.DOO/l28.DOO*(99.DOO*DSIN(6.DOO*DE)-12.000* 931 & DSIN(4.DOO*DE)-17.DOO*DSIN(2.DOO*0E)) 932 933 R(I)=M1*105.DOO/8.DOO*(33.DOO*DSIN(DE)**5-48.DOO* 934 & DSIN(0E)**3+16.0OO*DSIN(DE)) 935 GO TO 90 936 937 90 IF(DABS(P(I ) ).GT.PMAX)PDEG=DE 938 IF(DABS(P(I)).GT.PMAX)PMAX=DABS(P(I )) 939 IF(DABS(0(I ) ).GT.QMAX)QDEG=DE 940 IF(DABS(0(I)).GT.QMAX)QMAX=DABS(Q(I)) 941 IF(DABS(R(I)).GT.RMAX)RDEG=DE 942 IF(DABS(R(I)).GT.RMAX)RMAX=DABS(R( I)) 943 DE=DE+STEP 944 100 CONTINUE 945 RETURN 946 END 947 948 SUBROUTINE MODE(DE,AL,DXI,PH) 949 C 950 C Calculate the displacements at the corners of the Input unit 951 C area for various modes. At the moment only 1=8, m=8,7,6,5,4,3,2 952 C 1=6. m=6,5,4.3.2 , 1=7, m=7,6,5,4.3, are a v a i l a b l e 953 C 954 INTEGER L,M,LMODE,MMODE 955 REAL*8 DE,AL.DXI(3),PH,PI,BIDON(5),DCONST,M1 956 REAL*8 AMPL(3).PHIO(3).X,V 957 COMMON /E1/AMPL.PHI0/E4/BID0N,PI 95ff COMMON /E6/ LMODE,MMODE 959 M1=DFL0AT(MMODE) 960 M=IABS(MMODE) 961 IF((LMODE+1)-8)6,7.8 962 C L=8. M=+ OR - 8 963 8 IF(M+1-8)86.87.88 964 88 DXI( 1 ) = -AMPL(1)*OSIN(PHIO(1)+M1 *AL-PH*PI )* 965 & 2O27O25.DOO/128.DOO*(35.DOO+DC0S(8.0DOO*DE)-8.DOO* 966 & DC0S(6.DOO*DE)+28.0OO*DC0S(4.DOO*DE)-56.DOO* 967 5 DC0S(2.D0O*DE)) 968 DXI(2 ) = -AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI )* 969 & 2027025.DOO/128.D00*(-8.D00*DSIN(8.0OO*DE)+48.000* 970 & DSIN(6.DOO*DE)-112.DOO*DSIN(4.DOO*DE)+112.000* 971 a DSIN(2.DOO*DE)) 972 DXI(3) = -AMPL(3)*M1*DC0S(PHI0(3)+M1*AL-PH*P1 )* 973 a 2027025.DOO*(DSIN(DE)**7+6.DOO*DSIN(DE)) 974 RETURN 975 C L=8, M=+ OR -7 976 977 87 X=DCOS(DE) 978 Y=DSIN(DE) 979 DXI(1) = -AMPL(1)*DSIN(PHIO( 1)+M1 *AL-PH*PI)* 980 a 2027O25.0O0*(Y)*(-X**7+3.DOO*X**5-3.D0O*X**3+X) 981 982 DXI(2)=-AMPL(2)*DSIN(PHI0(2)+M1*AL-PH*PI)* 983 8 2027025.DOO*(Y**2*(7.D0O*X**6-15.D00*X**4+9.DOO* 984 , a X**2-1.DOO)+(-X**8+3.DOO*X**6-3.DOO*X**4+X**2)) 985 a 986 DXI (3) = -AMPL(3)*M1*DC0S(PHIO(3)+M1*AL-PH*PI)* I l l L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CCIOFRAS Page 18 987 & 2027025.DOO*(-X**7+3.DOO*X*»5-3.DOO*X**3+X) 988 RETURN 989 990 86 IF(M*1-6)84,85.886 991 992 C L=8, M= + OR - 6 993 994 886 X=DCOS(DE) 995 - Y=DSIN(DE) 996 DXI( 1 ) = -AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI)* 997 & 135135.DOO/2.DOO*(-15.DOO*X**8+46.DOO*X**6-48.DOO*X**4 998 & +18.DOO*X**2-1) 999 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 100O & 135135.DOO/2.000*Y*(120.DOO*X**7-276.DOO*X**5+192.DOO* 1001 & X**3-36.DOO*X) '1002 1003 0X1(3) =-AMPL(3)*M1*DCOS(PHI0(3)+M1 * AL-PH* PI)* 1004 & 135135.DOO/2.DOO*(14.DOO*Y**5-15.DOO*Y** 7) 1005 10O6 RETURN 1007 1008 C L=8. M=+ OR - 5 1009 1010 85 X=DCOS(DE) 101 1 Y=DSIN(DE) 1012 DXI(1)=-AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI)* 1013 & 135135.DOO/2.DOO*(Y)*(5.DOO*X**7-11.DOO*X**5+ 1014 & 7.DOO*X**3-X) 1015 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1016 & 135135.DOO/2.DOO*((5.DOO*X**8-11.000*X**6+7.DOO*X**4-1017 & X**2)-Y**2*(35.DOO*X**6-55.DOO*X**4+21.D00*X**2-1.DOO)) 1018 1019 DXI(3)=-AMPL(3)*M1*DC05(PHIO(3)+M1*AL-PH*PI)* 1020 & 135135.DOO/2.DOO*(5.DOO*X**7-11.DOO*X**5+ 1021 • & 7.D00*X**3-X) 1022 1023 RETURN 1024 1025 1026 84 IF (M+1-4)82,83,884 1027 1028 C L = 8, M=+ OR -4 1029 1030 884 X=DCOS(DE) 1031 Y=DSIN(DE) 1032 DXI(1)=-AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI)* 1033 & 10395.DOO/8.DOO*(65.000*X**8-156.DOO*X**6+118.DOO* 1034 & X**4-28.DOO*X**2+1 ) 1035 1036 DXI(2) = -AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI )* 1037 & - 10395*Y*(65.DOO*X**7-117.DOO*X**5 + 59.DOO*X**3-7.DOO*X) 1038 1039 DXI(3) = -AMPL(3)*M1*DCOS(PHIO(3) + M1 *AL-PH*PI)* 1040 & 10395.DOO/8.DOO*(65.D00*Y**7-104.DOO*Y**5+40.DOO*Y**3) 1041 1042 RETURN 1043 1044 C L = 8 i . M = + OR -3 112 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 19 1045 1046 83 X=DCOS(DE) 1047 Y=DSIN(DE) 1048 DXI( 1 ) = -AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI )* 1049 & 3465.DOO/8.DOO*Y*(-39.DOO*X**7+65.DOO*X**5-29.DOO* 1050 & X**3+3.DOO*X) 1051 1052 DXI(2) = -AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI )* 1053 & 3465.DOO/8.DOO*((-39 DOO*X**8+65.DOO*X**6-29.D00*X**4+ 1054 & 3D00*X**2)+(Y**2)*(273.DOO*X**6-325.DOO*X**4+87.DOO* 1055 & X**2-3.DOO)) 1056 DXI(3) = -AMPL(3)*M1*DC0S(PHIO(3)+M1 *AL-PH*PI)* 1057 & 3465.DOO/8.DOO*(-39.DOO*X* *7+65.DOO*X* *5-29.DOO* 1058 & X**3+3.DOO*X) 1059 1060 RETURN 1061 C L = 8 , M=+ OR -2 1062 1063 82 X=OCOS(DE) 1064 Y=DSIN(DE) 1065 DXI( 1 ) = -AMPL( 1 )*DSIN(PHIO(1)+M1*AL-PH*PI )* 1066 & 315.DOO/16.DOO*(-143.DOO*X**8+286.DOO*X**6-176.DOO* 1067 '& X**4+34.DOO*X**2-1.DOO) 1068 1069 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1070 & 315.DOO/4.DOO*(-Y)*(286.DOO*X**7+429.DOO*X**5-176.DOO* 1071 & X**3+17.DOO*X) 1072 1073 DXI(3) = -AMPL(3)*M1*DC0S(PHI0(3)+M1*AL-PH*PI )* 1074 & 315.DOO/16.DOO*(-143.DOO*Y**7+286.DOO*Y**5-176.DOO* 1075 & Y**3+32.DOO*Y) 1076 1077 RETURN 1078 C L=7, M = + OR - 7 1079 7 IF(M+1-7)75,76,77 1080 1081 77 X=DCOS(DE) 1082 Y=DSIN(DE) 1083 DXI(1)=-AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI)* 1084 & 135135.DOO*Y*(-X**6+3.DOO*X**4-3.DOO*X* *2+1 DOO) 1085 1086 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1087 & 135135.DOO*(-(Y**2)*(-6.DOO*X**5+12.DOO*X**3-6.DOO*X)+ 1088 & (-X**7+3.DOO*X**5-3.DOO*X**3+X)) 1089 1090 DXI(3)=-AMPL(3)*M1*DCOS(PHIO(3)+M1*AL-PH*PI)* 1091 & 135135.DOO*(-X**6+3.DOO*X**4-3.DOO*X**2+1.000) 1092 RETURN 1093 C L =7, M=+ OR -6 1094 1095 76 X=DCOS(DE) 1096 Y=DSIN(DE) 1097 DXI(1)=-AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI)* 1098 & 135135.DOO*(-X**7+3.DOO*X**5-3.DOO*X**3+X) 1099 1 100 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1 101 & 135135.DOO*(-Y)*(-7.DOO*X**6+15.DOO*X**4-9 D00*X**2+ 1 102 & 1.DOO) 113 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 20 1 103 1 104 DXI(3)=-AMPL(3)*M1*DC0S(PHI0(3)+M1*AL-PH*PI)* 1 105 & 135135.DOO*Y*(X**5-2.DOO*X**3+X) 1 106 RETURN 1 107 1 108 75 IF(M+1-5)73,74,775 1 109 1 1 10 C L=7, M=+ OR - 5 1111 1112 775 X=DCOS(DE) 1113 Y=DSIN(DE) 1114 DXI(1) = -AMPL(1)*DSIN(PHIO(1)+M1 *AL-PH*PI)* 1115 & 10395.DOO/2.DOO*(Y)*(13.DOO*X**6-27.DOO*X**4+15. DOO* 1 1 16 & X**2-1.DOO) 1117 DXI(2 ) = -AMPL(2)*DSIN(PHIO(2)+M1 *AL-PH*PI)* 1118 & 10395.DOO/2.DOO*(-(Y**2)*(78.DOO*X**5-108 D00*X**3+ 1 1 19 & 30.DOO*X)+(13.DOO*X**7-27.D00*X**5+15.DOO*X**3- X)) 1 120 1121 DXI(3)=-AMPL(3)*M1*DCOS(PHIO(3)+M1*AL-PH*PI)* 1 122 & 10395.DOO/2.DOO*(13.DOO*X**6-27.DOO*X**4+15.DOO* 1 123 & X**2-1.DOO) 1 124 1 125 RETURN 1 126 1 127 C L = ' 7, M=+ OR - 4 1 128 1 129 74 X=DCOS(DE) 1 130 Y=DSIN(DE) 1131 DXI(1)=-AMPL(1)*DSIN(PHIO(1)+M1*AL-PH*PI)* 1 132 & 3465.DOO/2.DOO*(13.DOO*X**7-29.D00*X**5+19.DOO* 1 133 & X**3-3.DOO*X) 1 134 1 135 DXI(2) = -AMPL(2)*DSIN(PHIO(2) + M1 *AL-PH*PI)* 1 136 & 3465.DOO/2.DOO*(-Y)*(91.DOO*X**6-145.DOO*X**4+57 . DOO* 1 137 & X**2-3.DOO) 1 138 1 139 DXI(3)=-AMPL(3)*M1*DC0S(PHI0(3)+M1*AL-PH*PI)* 1 140 & 3465.D0O/2.D0O*(Y)*(-13.D0O*X**5+16.DOO*X**3-3. DOO*X) 1 141 1 142 RETURN 1 143 1 144 1 145 73 IF (M+1-3)71,72,773 1 146 1 147 C L = 7 , , M=+ OR -3 1 148 1 149 773 X=DCOS(DE) 1 150 Y=DSIN(DE ) 1 151 DXI(1) = -AMPL(1)*DSIN(PHIO(1)+M1 *AL-PH*PI)* 1 152 & 315.DOO/8,DOO*(Y)*(-143.DOO*X**6+209.D00*X**4-69. DOO* 1 153 & X**2+3.DOO) 1 154 1 155 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1 156 & 315.DOO/8.DOO*(-(Y**2)*(-858.DOO*X**5+836.DOO*X**3-1 157 & 138.DOO*X) + (- 143.DOO*X**7+209.D00*X**5-69.DOO*X **3+ 1 158 & 3.DOO*X)) 1 159 1 160 DXI(3)=-AMPL(3)*M1*DCOS(PHIO(3)+M1*AL-PH*PI)* 114 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 21 1161 & 315.DOO/8.DOO*(-143.DOO*X**6+209.DOO*X**4-69.DOO* 1162 & X**2+3.DOO) 1163 1164 RETURN 1 165 1166 C L=7, M=+ OR -2 1 167 1168 72 X=DC0S(DE) 1169 Y =DSIN(DE) 1 170 DXI(1) = -AMPL(1)*DSIN(PHIO( 1)+M1 *AL-PH*PI ) 117 1 1172 DXI(2 ) = -AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI) 1 173 DXI(3)=-AMPL(3)*M1*DC0S(PHIO(3)+M1*AL-PH*PI) 1 174 1175 RETURN 1176 C L=7, M=+ OR -1 1 177 ' 1178 71 X=DCOS(DE) 1179 Y=DSIN(DE) 1180 DXI( 1 ) = -AMPL( 1 )*DSIN(PHIO(1)+M1*AL-PH*PI) 1181 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH»PI) 1182 DXI(3 ) = -AMPL(3 ) *M1*DCOS(PHIO(3)+M1*AL-PH*PI ) 1183 RETURN 1 184 1185 6 IF(M+1-6)664.65.66 1186 C L=6, M=+ OR -6 1 187 1188 66 DXI(1)=-AMPL(1)*DSIN(PHI0(1)+M1*AL-PH*PI)* 1189 & 11.DOO*945.DOO/32.DOO*(-DCOS(6.D00*DE)+6.DOO* 1190 & DC0S(4.DOO*DE)-15.DOO*DC0S(2.DOO*DE)+1O.DOO) 1191 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1192 & 66.000*945.DOO/32.DOO*(DSIN(6.DOO*DE)-4.DOO*DSIN( 4.DOO*DE) 1193 & +5.DOO*DSIN(2.DOO*DE)) 1 194 DXI(3) = -AMPL(3)*M1*DCOS(PHIO(3)+M1*AL-PH*PI)* 1195 & 11.DOO*945.DOO/16.DOO*(DSIN(5.DOO*DE)-5.DOO* 1196 & DSIN(3.DOO*DE)+10.DOO*DSIN(DE)) 1197 RETURN 1 198 C L = 6. M=+ OR -5 1199 65 DX 1(1) = -AMPL(1)*DSIN(PHIO(1)+M1 *AL-PH*PI ) * 1200 a 11.000*945.DOO/32.DOO*(DSIN(6,DOO*DE)+4,D00*DSIN(4.DOO*DE) 1201 a +5.DOO*DSIN(2.DOO*DE)) 1202 DXI(2) = -AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI )* 1203 a 11.DOO*945.DOO/32.DOO*(6.DOO*DCOS(6.DOO*DE)+16.DOO 1204 a *DC0S(4.D0O*DE)+10.D00*DC0S(2.D0O*DE)) 1205 DXI(3)=-AMPL(3)*M1*DC0S(PHI0(3)+M1*AL-PH*PI)* 1206 a 11.D0O*945.DOO/16.DOO*(DC0S(5.DOO*DE)-3.DOO* 1207 a DC0S(3.DOO*DE)+2.DOO*DC0S(DE)) 1208 RETURN 1209 664 IF(M+1-4)62,63.64 1210 1211 C L=6,M=+ OR -2 1212 62 DXI(1 ) = -AMPL(1 )*DSIN(PHIO(1)+M1*AL-PH*PI)* 1213 8 105.DOO/256.DOO*(10.DOO+17.DOO*DCOS(2.DOO*DE)+6.DOO* 1214 a DC0S(4.DOO*DE)-33.DOO*DC0S(6.DOO*DE)) 1215 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1216 a 105.DOO/128.DOO*(99.DOO*DSIN(6.DOO*DE)-12.DOO* 1217 a DSIN(4.DOO*DE)-17.DOO*DSIN(2.DOO*DE)) 1218 DXI(3)=-AMPL(3)*M1*DCOS(PHIO(3)+M1*AL-PH*PI)* 115 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CCid=FRAS Page 22 1219 & 105.DOO/8.DOO*(33.DOO*DSIN(DE) **5-48.DOO*DSIN(DE)**3 1220 & +16.D00*DSIN(DE)) 1221 RETURN 1222 1223 1224 C L=G.M=+ OR -3 1225 1226 63 DXI(1)=-AMPL(1)*DSIN(PHI0(1)+M1*AL-PH*PI)* 1227 & 315.DOO/64.D00*(- 11.DOO*DSIN(6.D00*DE)+12.DOO* 1228 & DSIN(4.DOO*DE)+9.DOO*DSIN(2.DOO*DE)) 1229 1230 DXI(2)=-AMPL(2)*DSIN(PHIO(2)+M1*AL-PH*PI)* 1231 & 1890.DOO/64.DOO*(-11.DOO*DC0S(6.DOO*DE)+8.DOO* 1232 & DCOS(4.DOO*DE)+3.DOO*DC0S(2.DOO*DE)) 1233 1234 DXI(3)=-AMPL(3)*M1*DC0S(PHIO(3)+M1*AL-PH*PI)* 1235 & 315.DOO/32.DOO*(55.DOO*DCOS(DE)+DCOS(3.DOO*DE)-11.DOO* 1236 & DC0S(5.DOO*DE)) 1237 1238 RETURN 1239 1240 C L=6, M=+ OR -4 1241 64 DXI(1 ) = -AMPL( 1 )*DSIN(PHIO(1)+M1 *AL-PH*PI )* 1242 & 9.DOO*105.DOO/64.DOO*(11.DOO*DC0S(6.DOO*DE)-26.DOO 1243 & *DC0S(4.0O0*DE)+5.DO0*DC0S(2.DOO*DE)+1O.DOO) 1244 DXI(2) = -AMPL(2)*DSIN(PHI0(2)+M1*AL-PH*PI )* 1245 & 9.DOO*105.DOO/64.DOO*(-66.D00*DSIN(6.DOO*DE)+104.DOO 1246 & *DSIN(4.DOO*DE)-10.DOO*DSIN(2.DOO*DE)) 1247 DXI(3) = -AMPL(3)*M1*DCOS(PHIO(3)+M1 *AL-PH*PI)* 1248 & 9.DOO*105.DOO/32.DOO*(-11.0OO*DSIN(5.DOO*DE)+15.DOO* 1249 & DSIN(3.DOO*DE)+10.DOO*DSIN(DE)) 1250 RETURN 1251 END 1252 1253 SUBROUTINE PARAM(VRMAX.XRMAX,K,PMAX,OMAX,RMAX,AMPL,SIG,SIGTD, 1254 &VTIDE,TIDE ) 1255 C This Subroutine calculates the amplitudes and angular frequenc 1256 C the mode, and v e l o c i t y and displacement parameters being used 1257 REAL*8 VRMAX,XRMAX,K,PMAX,OMAX,RMAX.SIGTD 1258 REAL*8 AMPL(3),SIG(3),TIDE(3) 1259 XRMAX=-VRMAX/SIG(1) 1260 AMPL(1)=-XRMAX/PMAX 1261 AMPL(2)=(K*VRMAX)/(SIG(1)*QMAX) 1262 AMPL(3)=(K*VRMAX)/(SIG(1)*RMAX) 1263 SIG(2)=SIG(1) 1264 SIG(3)=SIG(1) 1265 TIDE(1)=VTIDE/(3.DOO*SIGTD) 1266 TIDE(2)=(K*VTIDE)/(SIGTD*6.DOO) 1267 TIDE(3) = (K*VTIDE )/(SIGTD*6.DOO) 1268 RETURN 1269 END 1270 1271 C 1272 SUBROUTINE POLCA(MP.MC,J) 1273 C 1274 C Convert from polar to cartesian coordinates 1275 C 1276 REAL*8 MC(3,J).MP(3.d) 116 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 23 1277 DO 1 K=1.J 1278 MC(3,K)=MP( 1 ,K)*DC0S(MP(2,K)) 1279 MC(2,K)=MP(1,K)*DSIN(MP(2.K)) 1280 MC(1,K)=MC(2,K)*DC0S(MP(3,K)) 1281 1 MC(2.K)=MC(2.K)*DSIN(MP(3.K)) 1282 RETURN 1283 END 1284 C 1285 SUBROUTINE RADIUS(DEC.RAD) 1286 C ,1287 REAL*8 DEC,RAD,0BL,RP0L,BID0N(2) 1288 COMMON /E4/BID0N,OBL,RPOL 1289 RAD=RP0L/DS0RT(1.D+00+OBL*(OBL-2.D+OO)*DSIN(DEC)*DSIN(DEC)) 1290 RETURN 1291 END 1292 1293 C 1294 SUBROUTINE ROLL(MAT) 1295 C 1296 REAL*8 MAT(3.4),MEM(3) 1297 DO 1 1=1.3 1298 1 MEM(I)=MAT(1,1) 1299 DO 2 J=1.3 1300 DO 2 1=1.3 1301 2 MAT(I,J)=MAT(I,J+1) 1302 DO 3 1=1,3 1303 3 MAT(I,4)=MEM(I) 1304 RETURN 1305 END 1306 SUBROUTINE SMEAR(Y.C0NV2,NN) 1307 C Program reads In a model l i n e , generates an appropriate smearing 1308 C gauss 1 an, and convolves the two using Subroutine APCONV 1309 C APCONV courtesy of Gerry Grieves 1310 C 131 1 INTEGER M2.NN 1312 REAL*8 LAMO,DLAM,NSTEP 1313 REALM SIGMA. AREA 1 . ARE A2 , N. SDLAM, C 1314 REAL*4 DATA(512).GAUSS(512).C0NV(512).C0NV2(512),Y(512) 1315 COMMON /E3/ENDS.LAMO.NSTEP 1316 COMMON /E7/DLAM 1317 SDLAM=SNGL(DLAM) 1318 C 1319 C R e c t i f y the input l i n e p r o f i l e 1320 C 1321 DO 10 I=1,NN 1322 DATA(I) = 1.OOEOO-Y( I ) 1323 AREA1=AREA1+DATA(I) 1324 10 CONTINUE 1325 1326 c 1327 c Generate a Gaussian of FWHM=.3 angstroms, cut off occurs when the 1328 c value is less than 1.00d-06 1329 c 1330 SIGMA=7.2135E-03 ! 1.27398E-02 1331 C=O.OOEOO 1332 DO 20 1=1.512 1333 GAUSS(I)=EXP(-C**2/(2.E00*SIGMA**2 ) ) 1334 C=C+SDLAM 117 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CC1d=FRAS Page 24 1335 IF(GAUSS(I).LT.1.OOE-04)G0 TO 25 1336 20 CONTINUE 1337 25 X=FLOAT(I) !X=0of points In the gaussian p r o f i l e 1338 C=O.OOEOO 1339 DO 30 d=1,100 1340 GAUSS(d)=EXP(-(C-(X*SDLAM))**2/(2*SIGMA**2)) 134 1 C=C+SDLAM 1342 IF(GAUSS(J).LT.1.OOE-04)GAUSS(J)=0.OOEOO 1343 30 CONTINUE 1344 1345 C 1346 C Ca l l APCONV to c a l c u l a t e the convolution 1347 C 1348 M=(2*I) 1349 CALL APCONV(DATA,512,GAUSS.M.C0NV.M2) 1350 1351 C 1352 C 1353 DO 55 d=1,M2 1354 AREA2 = AREA2+C0NV(d) 1355 55 CONTINUE 1356 A=AREA1/AREA2 1357 DO 60 1=1,M2 1358 C0NV2(I)=C0NV(I)*A 1359 C0NV2(I)=1.0OEOO-C0NV2(I) 1360 60 CONTINUE 1361 RETURN 1362 END 1363 1364 1365 1366 1367 1368 1369 1370 1371 C 1372 ' SUBROUTINE SPEC(Z.RS) 1373 C 1374 C Move the i n t r i n s i c l i n e p r o f i l e according to the c a l c u l a t e d v e l o c i t y 1375 C and add i t to the tot a l l i n e p r o f i l e . Note that l i n e a r i n t e r p o l a t i o n 1376 C is used between the wavelength bins. 1377 C 1378 INTEGER NSTEP,ENDS(2),MAX,MIN 1379 REAL*8 Z,RS,LAMO.DLAM,KA,B 1380 REAL*8 SPECTO(SOO).SPECTF(500),LAM(500),C 1381 COMMON /E3/ENDS.LAMO.NSTEP 1382 COMMON /E5/SPECT0,SPECTF.LAM 1383 COMMON /E7/DLAM 1384 1385 C=2.99792458D+08 1386 B=Z*LAMO/(C*DLAM) 1387 . KA=1.D+OO+Z/C 1388 IF (Z) 2,3.4 1389 3 I1=ENDS(1) 1390 I2=ENDS(2) 1391 DO 5 1 = 1 1 , 1 2 1392 5 SPECTF(1 + 1)=SPECTF(1 + 1 ) + RS*SPECTO(1+1) 118 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8. 1985 for CCid=FRAS Page 25 1393 GO TO 6 1394 2 MAX=IDINT(B+KA*DBLE(FLOAT(NSTEP))) 1395 IF (MAX.LT.ENDS(2)) ENDS(2)=MAX 1396 GO TO 8 1397 4 MIN=1+IDINT(B) 1398 IF (MIN.GT.ENDS(1)) ENDS(1)=MIN 1399 8 I1=ENDS(1) 1400 I2=ENDS(2) 1401 NN=I2-I1+1 1402 DO 10 1=1,NN 1403 d=IDINT(1.D+00+(DBLE(FLOAT(I+I1))-1 D+00-B)/KA) 1404 10 SPECTF(I + I1) = SPECTF(I + I1)+RS*(SPECT0(d)*LAM( d+1) 1405 Q-SPECT0(d+1)*LAM(d)+LAM(I+I1)*(SPECTO(J+1) 1406 0-SPECTO(J))/KA)/DLAM 1407 6 CONTINUE 1408 RETURN 1409 END 14 10 C 14 11 C 1412 SUBROUTINE SPEED(VP.VXI,VXIP) 1413 C 1414 C Calculate the v e l o c i t y of the unit area due to the nonradial 1415 C o s c i l l a t i o n s 1416 C 1417 INTEGER IL 1418 REAL*8 VP(3).VXI(3),VXIP(3).OME,MU.DXI(3) 1419 REAL*8 SIG(3).0X11(3).SIGTD 1420 1421 COMMON /E4/MU.0ME/E2/SIG 1422 COMMON /TIDE/SIGTD,TIDE,TPHS 1423 1424 DO 1 1=1.3 1425 1 VXIP(I)=0.D+OO 1426 1427 CALL MODE(VP(2),VP(3).DXI,O.5D+00) ! vel contribution from o s c i l l a t i o n 1428 DO 2 1=1,3 1429 2 VXIP(I ) = VXIP(I) + SIG(I)*DXI(I) 1430 1431 CALL ATI0E(VP(2),VP(3),DXI1,0.5D+00) ! vel cont from the t i d a l e f f e c t 1432 DO 3 1=1.3 1433 3 VXIP(I)=VXIP(I)+SIGTD*DXI1(I) 1434 1435 VXIP(3) = VXIP(3)+0ME*VP(1)*DSIN(VP( 2)) 1436 VXI(1)=DC0S(VP(3))*(VXIP(1)*DSIN(VP(2))+VXIP(2)* 1437 0DC0S(VP(2)))-DSIN(VP(3))*VXIP(3) 1438 VXI(2)=DSIN(VP(3))*(VXIP(1)*DSIN(VP(2))+VXIP( 2)* 1439 0DC0S(VP(2)))+DCOS(VP( 3))*VXIP(3) 1440 VXI(3)=DC0S(VP(2))*VXIP(1)-DSIN(VP(2))*VXIP(2) 1441 RETURN 1442 END 1443 C 1444 SUBROUTINE SURF(MAT,S,VS.VC) 1445 C 1446 C Calculate displacement of the unit surface due to the NR o s c i l l a t i o n 1447 C 1448 REAL*8 MAT(3.4).S.VS(3),VC(3).VA(3),VB(3),VSB(3),M1 1449 RE AL*8 M2 , RM 1450 DO 1 1=1,3 119 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CC1d=FRAS Page 26 1451 VC(I)=0.D+00 1452 1 VS(I)=0.D+00 1453 DO 2 1=1,4 1454 DO 3 d=1,3 1455 VA(d)=MAT(J,3)-MAT ( J . 2) 1456 3 VB(J)=MAT(J,1)-MAT(J,2) 1457 CALL CR0SSC(VA,VB,VSB) 1458 DO 4 J=1.3 1459 4 VS(J)=VS(J)+VSB(J)/4.D+00 1460 2 CALL ROLL(MAT) 1461 S=O.D+00 1462 DO 5 1=1.3 1463 5 S=S+VS(I)*VS(I) 1464 S=DSQRT(S) 1465 DO 6 1=1,3 1466 6 VS(I)=VS(I)/S 1467 DO 7 1=1,4 1468 DO 7 J=1,3 1469 7 VC(J)=VC(J)+MAT(d,I)/4.D+OO 1470 M1=O.D+00 1471 DO 8 1=1,4 1472 DO 8 J=1,3 1473 8 M1=M1+MAT(J,I)*MAT(J,I)/4 D+OO 1474 M2=0.D+00 1475 DO 9 1=1.3 1476 9 M2=M2+VC(I)*VC(I) 1477 RM=DS0RT(M1/M2) 1478 DO 10 1=1,3 1479 10 VC(I)=VC(I)*RM 1480 RETURN 1481 END 1482 1483 1484 C > 1485 SUBROUTINE TRACE(NN) 1486 C 1487 C Outputs the f i n a l l i n e p r o f i l e . Also c a l l s SMEAR and BINDAT to 1488 C produce a lower r e s o l u t i o n l i n e p r o f i l e . 1489 C 1490 1491 INTEGER ENDS(2).NN,NINDEX 1492 REAL*8 SPECTO(50O).SPECTF(500) 1493 REAL*8 LAM(50O) 1494 REAL*4 X(512 ) ,Y(5 12 ) .CONV(512),MAX 1495 REAL*4 BIN(512) 1496 COMMON /E3/ENDS,LAMO,NSTEP 1497 COMMON /E5/SPECT0.SPECTF,LAM 1498 COMMON /E7/DLAM 1499 1500 100 FORMATC1 SYNTHETIC SPECTRUM:'/(' LAMBDA= '.F15.5. 1501 0' NM',10X,'INTENSITY' '.E18.7)) 1502 1503 I1=ENDS(1) 1504 I2=ENDS(2) 1505 MAX=0.0 1506 1507 WRITE ( 16.100)(LAM(1 + 1).SPECTFC1 + 1) . 1 = 1 1 . 12) 1508 C 120 L i s t i n g of PROGRAM.FOR at 18:40:14 on MAR 8, 1985 for CCi d= FRAS Page 27 1509 C Determine the maxium in the spectrum and then normalise i t 1510 C 1511 DO 1 1=11.12 1512 IF (MAX.LT.SNGL(SPECTF(I+1))) MAX=SNGL(SPECTF(I+1)) 1513 1 CONTINUE 1514 DO 2 I=1.NN 1515 Y(I) = SNGL(SPECTF(I 1 + 1))/MAX 1516 2 X(I)=SNGL(LAM(I 1 + 1)) 1517 1518 1519 1520 1521 C 1522 C C a l l SMEAR-Smear w i l l output a f i l e containing the smeared 1523 C l i n e , and one containing the o r i g i n a l model l i n e 1524 C Y contains the model line.CONV contains the smeared l i n e 1525 C 1526 1527 NINDEX=300 1528 N11=80 1529 1530 CALL SMEAR(Y.CONV.NN) 1531 CALL BINDAT(CONV.BIN.NINDEX) 1532 1533 BIN( 1 )=.90 1534 WRITE(19,500) N11 1535 WRITE(19.501)(BIN(I).1=1,N11) 1536 WRITE(6.*)(BIN(I).I=1.N11) 1537 Y(1)=.90 1538 WRITE(17.500) NINDEX 1539 WRITE(17,501)(Y(I),1=1.NINDEX) 1540 1541 CONV(1)=.90 1542 WRITE(18,500) NINDEX 1543 WRITE(18,5O1)(C0NV(I),I=1.NINDEX) 1544 1545 500 FORMAT(15) 1546 501 F0RMAT(5F10.5) 1547 RETURN 1548 END 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085599/manifest

Comment

Related Items