Парадоксы при аналоговом и математическом моделировании параметрического резонанса

Парадоксы при аналоговом и математическом моделировании параметрического резонанса

Аннотация

Работа является продолжением экспериментальной работы, которая докладывалась на конференции в МФТИ в прошлом году. На этот раз основное внимание уделено интерпретации результатов эксперимента. Показано, как с помощью простого механического осциллятора можно экспериментально изучать процессы в колебательном контуре, которые сопровождаются непосредственным превращением механической энергии в энергию электромагнитного поля. Результаты экспериментов подтверждаются численным решением дифференциального уравнения, которое описывает наблюдаемое физическое явление.

В качестве простейшего устройства для аналогового моделирования выбран маятник на трифилярном подвесе. Показано, что физические процессы в колебательном контуре и в выбранном механическом осцилляторе описываются одинаковыми дифференциальными уравнениями,

Особое внимание уделяется силе Кориолиса. Показано, что параметрическое изменение амплитуды колебаний в маятнике с тремя подвесами обусловлено силой Кориолиса, действующей на движущиеся дополнительные грузы в процессе колебаний.

Сопоставляются уравнения для колебательного контура с изменяющейся индуктивностью и уравнение для механического осциллятора. Показано, что при параметрическом изменении амплитуды колебаний момент инерции в маятнике на трех подвесах и индуктивность катушки в колебательном контуре играют одну и ту же роль.

Полученные теоретические результаты подтверждаются экспериментом и  численными расчетами. Работа может быть использована при создании новых лабораторных работ по физике.

1. Введение

  Простейшим устройством для превращения механической энергии непосредственно в электромагнитную энергию является колебательный контур. Физическое  явление, которое при этом используется, называется параметрическим резонансом. В отличие от обычного резонанса, обусловленного действием периодической вынуждающей силы или периодического возбуждающего источника, например, источника ЭДС, подпитывающего колебательный контур, параметрический резонанс наблюдается при изменении параметров осциллирующей системы [1].  Такими изменяющимися параметрами в колебательном контуре являются емкость и индуктивность.

Интересно, что каждый человек использует параметрический резонанс уже в раннем детстве, когда раскачивает качели без помощи извне, приседая и вставая на них в определенные моменты времени. 

Для наблюдения параметрического резонанса в колебательном контуре требуются специальные  приспособления, позволяющие быстро выдвигать сердечник из катушки или диэлектрическую пластину из конденсатора и вдвигать их обратно в определенные моменты времени. Такая установка, в которой механическая энергия, затрачиваемая на вытаскивание ферромагнитного сердечника из катушки, превращается в электрическую энергию,  была создана в двадцатых годах прошлого века  и  [1].

Оказывается, изучение параметрического резонанса в колебательном контуре возможно без самого контура, если воспользоваться  простым  механическим  осциллятором. Имеется в виду реальный физический эксперимент, а не математическое моделирование, которое, казалось бы,  отодвинуло на второй план аналоговое физическое моделирование с помощью лабораторных моделей, свойства которых описываются уравнениями, аналогичными уравнениям изучаемого объекта. Обладая рядом преимуществ, аналоговое моделирование позволило получить результаты принципиальной важности, например, при исследовании волн в плазме  и в интерпретации вихревых структур в атмосферах планет. Медленно дрейфующее вокруг оси планеты Большое Красное Пятно Юпитера (БКПЮ)  наблюдали чуть ли не в  кастрюле с водой [2].

В настоящей работе используется своя “кухня” – моделирование с помощью “гвоздя и веревки”. В качестве простейшего устройства для аналогового моделирования выбран маятник на трифилярном подвесе. Как будет показано ниже, физические процессы в колебательном контуре и в выбранном механическом осцилляторе описываются одинаковыми дифференциальными уравнениями,

  2. Вывод уравнения свободных колебаний маятника

Что представляет собой маятник на трифилярном подвесе?  Это три лежащих в горизонтальной плоскости одинаковых стержня (три “гвоздя”) длиной d, соединенных своими  концами в одной точке О и подвешенных на трех нитях (на трех “веревках”), привязанных к стержням на  расстоянии  b от точки О. Угол между двумя соседними стержнями составляет 120о, длина каждой нити равна l.  масса  каждого стержня равна M/3. На каждом стержне  находится  дополнительный груз, массой m/3. Грузы синхронно  перемещаться по стержням, не приближаясь к точке О на расстояние, меньшее, чем b. Расстояние от точки  О до грузов обозначим x(t) (x(t)>b). Угол поворота стержней вокруг вертикали ОО1 обозначать . Наряду с геометрическими параметрами l, d, b и x  отношение масс  M/m и ускорение свободного падения g также  можно рассматривать как  параметры маятника.

Если маятник совершает свободные колебания при неподвижных  дополнительных  грузах, его энергия Е остается постоянной, поэтому уравнение колебаний можно получить, приравнивая нулю производную от Е по времени. Энергия маятника Е=Ек+Ер - это сумма кинетической энергии Ек вращательного движения стержней вокруг оси ОО1 и потенциальной энергии Ер  стержней, связанной с их  перемещением во время колебаний по вертикали вдоль той же оси ОО1 на величину .

.  (1).

Здесь кинетическая энергия выражена через момент инерции стержней и дополнительных грузов, - производная по времени от угла отклонения маятника .

Для исключения  из уравнения (1) решением геометрической задачи легко находится связь между углами  и  :

.  (2)

После исключения  и дифференцирования получим

,  (3)

Пренебрегая малой величиной  ,  получим уравнение колебаний

  (4)

  с циклической частотой

    (5) 

При малых углах   в равнении (4) можно положить .

А каким уравнением будут описываться колебания маятника, если дополнительные грузы придут в движение? Прежде, чем ответить на этот вопрос, сформулируем задачу для маятника на трифилярном подвесе, которая является аналогом для  задачи о колебательном контуре с резко изменяющейся индуктивностью. 

  3. Механический аналог задачи о колебательном контуре

Задача 1. Описанный выше маятник на трех подвесах совершает незатухающие колебания малой амплитуды около вертикальной оси.  В момент, когда стержни маятника проходят положение равновесия, дополнительные грузы одновременно быстро перемещаются с одинаковой скоростью от концов стержней к точкам подвеса. В состоянии наибольшего отклонения стержней  от положения равновесия дополнительные грузы быстро возвращаются обратно. Найти отношение увеличения угловой амплитуды маятника к первоначальной амплитуде после двукратного выполнения описанной процедуры.

Как нетрудно будет заметить, задача 1 чуть ли не дословно повторяет следующую задачу о параметрическом резонансе в колебательном контуре, которая предлагалась на одной из олимпиад  МФТИ. Методы решения таких задач обсуждаются в статье [3].

Задача 2. В колебательном контуре с индуктивностью L и емкостью C происходят незатухающие колебания. В момент, когда заряд конденсатора достигает максимального значения,  в катушку индуктивности быстро вставляют ферромагнитный сердечник, что повышает ее индуктивность до значения . В момент достижения током максимального значения стержень быстро выдергивают. Эту процедуру,  приводящую к увеличению амплитуды колебаний заряда,  периодически повторяют.  Найти отношение  увеличения амплитуды заряда конденсатора  при двукратном выполнении описанной процедуры к первоначальному заряду конденсатора

  Решение задачи 2 легко получить с помощью закона сохранения энергии, если предположить, что магнитный поток, пронизывающий катушку не успевает измениться за время вытаскивания сердечника из катушки. Для того чтобы сердечник вытащить из катушки, требуется преодолеть силу, с которой магнитное поле,  втягивает в катушку этот сердечник, На изменение индуктивности при вытаскивании сердечника из катушки  затрачивается  механическая  работа, за счет которой с увеличением  амплитуды тока  происходит  прирост  электромагнитной  энергии в  контуре.

Пусть I1 – ток в катушке перед выдергиванием сердечника, а I2 – ток в катушке сразу после того, когда сердечник убрали. Из сохранения магнитного потока  µLI1 =LI2  следует  увеличение тока в µ раз, I2=µI1.

Теперь легко убедиться, что описанная в задаче манипуляция с сердечником  приводит к увеличению энергии катушки тоже в µ раз,

  (6) 

Очевидно, повторение описанной процедуры два раза приведет к увеличению максимальной энергии  катушки, а значит, и энергии контура в µ2 раз. Если Q – первоначальная амплитуда заряда конденсатора и ∆Q – увеличение амплитуды заряда после увеличения энергии контура в µ2 раз, то из отношения энергий конденсатора

  следует : .

Если ввести величину увеличения индуктивности катушки  ∆L=µL-L, то последняя формула принимает симметричный вид.

  (7)

Отношение увеличения амплитуды заряда конденсатора к его  первоначальному заряду  равно отношению изменения индуктивности катушки  к меньшему значению индуктивности (индуктивности катушки без сердечника).

Для сравнения приведем решение задачи 1.  Из закона сохранения момента импульса: найдем  отношение угловых скоростей  , где

, .

При однократном выполнении описанной в задаче  процедуры полная энергия маятника возрастет в раз,  а при двукратном в раз.

, откуда .

Наконец, отношение изменения амплитуды к первоначальной амплитуде равно отношению изменения момента инерции маятника к меньшему  значению момента инерции

  (8)

Полученный результат такой же, что и в задаче о колебательном контуре, При параметрическом изменении амплитуды колебаний момент инерции в маятнике на трех подвесах и индуктивность катушки в колебательном контуре играют одну и ту же роль.

Появляется соблазн рассмотреть аналогичную задачу.

Задача 3. В колебательном контуре с индуктивностью L и емкостью C происходят незатухающие колебания. В момент, когда ток в катушке  достигает максимального значения,  в конденсатор быстро вставляют диэлектрическую пластину с диэлектрической проницаемостью  ε, что повышает его емкость до значения εС. В момент достижения зарядом конденсатора максимального значения пластину  быстро выдергивают.  Эту процедуру,  приводящую к увеличению амплитуды колебаний заряда,  периодически повторяют.  Найти отношение  увеличения амплитуды заряда конденсатора  при двукратном выполнении описанной процедуры к первоначальному заряду конденсатора.

Для этой задачи тоже можно найти механический аналог в маятнике на трех подвесах, если использовать только один дополнительный груз, который подвешивается в точке О.

4. Сила Кориолиса в маятнике на трех подвесах с подвижными грузами

Сначала получим уравнение для колебательного контура с изменяющейся индуктивностью L(t) и с постоянной  емкостью C. . Для этого приравняем напряжение на конденсаторе и ЭДС индукции в катушке,

  (9) 

Как только что было показано, момент инерции в маятнике на трех подвесах является механическим аналогом индуктивности . Отправляясь от уравнения (4), это позволяет записать уравнение колебаний маятника с подвижными грузами в виде

.  (10)

Здесь- механический аналог емкости:

.  (11)

Вычислив производную,  уравнение (10) можно записать в виде

  (12)

Второй член левой части уравнения (12) представляет собой момент силы Кориолиса . Это указывает на то, что параметрическое изменение амплитуды колебаний в маятнике с тремя подвесами обусловлено силой Кориолиса, действующей на движущиеся грузы.

Окончательно уравнение (12) запишем в виде

.  (13)

Уравнения (10) и (13) дополняется начальными условиями, . Обозначив угловую скорость , получим систему  двух дифференциальных уравнений первого порядка, которая решается численно.

  5. Опыты  с маятником

Смогут  ли три муравья одинаковой массы, двигаясь по стержням маятника, изменить его амплитуду  колебаний без помощи извне?  Конечно, послушные муравьи – это мечта экспериментатора. В нашем маятнике вместо них используются три дополнительных груза одинаковой массы, которые могут скользить по стержням. К каждому грузу привязывается  тонкий легкий резиновый жгут. Растягивая жгут, грузы смещают на край стержней и удерживают в таком положении с помощью тонкой нити, связанной с концами всех трех стержней. В нужный момент нить пережигается, и грузы под действием силы упругости одновременно приходят в движение. Если же нить перерезать ножницами, то избежать воздействия извне  не удается.

Двигаясь по стержню, муравьи совершают работу по преодолению центробежной силы инерции. По закону сохранения энергии величина этой работы равна приросту  энергии маятника. Для того, чтобы прирост энергии был максимальным, возвращаться назад  муравьи должны  тогда, когда  угловая скорость маятника минимальна и центробежная сила инерции мала.

Далее  результаты опытов подтверждаются численными расчетами.

6. Результаты численных расчетов

Далее  результаты опытов подтверждаются численными расчетами. Расчеты проводились для  l=0.7, d=0.3, b=0.1,  0.1<x<0.3 и  M/m=0.2.  Можно говорить лишь об условном  периоде колебаний Т=2.8 с. Колебания ведь не являются периодическими. Скорость  движения грузов изменялась так, что производная от момента инерции по времени оставалась постоянной в течение всего времени движения. Для этого указывалось время в долях периода, за которое грузы проходят расстояние (d – b). Методом прогноза и коррекций [4] решалась система из двух  уравнений, заменяющая одно уравнение (13).

с начальными условиями , . Во всех вариантах расчета полагалось . При счете контролировалось выполнение закона сохранения энергии.

Рис. 1.  Увеличение угла отклонения маятника в радианах со временем.  , , .

Рис. 2.  Зависимости угловой скорости в с-1 от времени в секундах.  , , .

Рис. 3  Угол в радианах, угловая скорость в  с-1. Фазовый портрет маятника для числовых данных, которые были использованы в рис. 1. и рис. 2.  , , .

Полученные численные результаты, представленные на рисунках 1 – 3, совпадают с аналитическими результатами, полученными выше при решении задачи 1 в предположении, что  грузы проходят расстояние (d – b) быстро, а амплитуда колебаний невелика.

На рис. 2 показано резкое увеличение амплитуды угловой скорости в моменты времени, когда стержни проходят положение равновесия. На рис. 3 на фазовом портрете хорошо просматриваются прямолинейные участки, соответствующие этим временам движению дополнительных грузов.

На следующих трех рисунках  представлены результаты расчетов, которые отличаются от результатов предыдущих трех рисунков только большим начальным углом отклонения. Колебания прекращаются, маятник начинает закручиваться в одну сторону.

  Рис. 4.  Зависимость угла отклонения маятника в радианах от времени в секундах для , ,

.Как показано на рис. 4,  угол отклонения вдруг начал увеличиваться. Это означает, что маятник стал вращаться в одну сторону, т. е.. колебания прекратились. Это явление также наблюдалось в опытах с маятником на двух бифилярных подвесах.

Рис. 5.  Зависимости угловой скорости в с-1 от времени в секундах для , , ..На рис. 5 показано, что угловая скорость. маятника изменяется по величине, но не меняет направления.  Но этого в экспериментах не наблюдается. Оказывается, наша модель описывает в этом случае не маятник, а что - то похожее на велосипедное колесо с золотником.

Рис. 6  Фазовый портрет маятника для числовых данных, которые были использованы в рис. 4. и рис. 5,  для , , . Угол в радианах, угловая скорость в  с-1.

Рис. 7  Фазовый портрет маятника для числовых данных, которые были использованы в рис. 1 и рис. 2,  но для скорости движения грузов в 4 раза меньшей,  , , .. Угол в радианах, угловая скорость в  с-1.

Как показано на рис. 7, с увеличением времени движения грузов резких изменений на фазовом портрете, таких, как на рис. 3, не наблюдается.

Выводы

1. Для изучения параметрического резонанса в колебательном контуре  предлагается проводить реальные физические эксперименты с его механическим аналогом – маятником на трех подвесах с подвижными дополнительными грузами на его стержнях.

2. Показано, что не существует принципиальной разницы между параметрическими резонансами в колебательном контуре и в маятнике с трифилярным подвесом.

3.  Показано, что параметрическое изменение амплитуды колебаний в маятнике с тремя подвесами обусловлено силой Кориолиса, действующей на движущиеся дополнительные грузы в процессе колебаний.

4. Полученные теоретические результаты подтверждаются экспериментом и  численными расчетами.

5. Работа может быть использована при создании новых лабораторных работ по физике.

В заключение остается выразить благодарность научному руководителю

Литература

1.   Эволюция понятия резонанса. Успехи физических наук. Т. ХХХI, вып. 4,  1947 г,  с. 447.

2. .   Вихри Россби и спиральные структуры: Астрофизика и физика плазмы в опытах на мелкой воде. –М.: Наука. 1990. -240 с.

3. В.  Можаева  Катушки индуктивности в электрических цепях  ”Квант” № 4 , 1998 г.

4. Мак- , исленные метолы и программирование на фортране.  Мир 1969.

Программа на фортране

        PROGRAM IGOR

  COMMON/CCC/KTURN, KSW, DD, ALP, ALM, TIME, BB, CC, CD3,AmM, XNP1,T0,XX, ALS

  6 FORMAT(1X,5F10.4,2I7,E12.4)

  7 FORMAT(1X,9F10.6)

cccccccccccccccccccccccccccccccccccccccccccccccc

  KTURN=0

        KSW=1

  DD=0.30

        BB=0.1

        AmM=0.1

        CD3=DD*DD/3.

cccccccccccccccccccccccccccccc

        Am=AmM*0.2

cccccccccccccccccccccccccccccc

        AL=0.7

        GG=10.0

        GL=GG/AL

        B2M=BB*BB*(1.0+AmM)

        WW=SQRT(B2M/(CD3+AmM*BB*BB)*GG/AL)

        PI=3.1415926

        TT=2.*PI/WW

        T0=TT/20.0        *4.0

cccccccccccccccccccccccccccccccccc

        VM=Am*(DD-BB)/T0

                VV=(DD-BB)/T0

cccccccccccccccccccccccccccccccccc

        CC=AL/(BB**2*GG*(1.0+AmM))

       ALS=AmM*(DD**2-BB**2)/T0

       ALP=CD3+AmM*DD*DD

               ALM=CD3+AmM*BB*BB

       OTBET=AmM*(DD*DD-BB*BB)/(1./3.*DD*DD+AmM*BB*BB)

        WRITE (*,7) TT, T0,OTBET, WW, CC, ALS 

cccccccccccccccccccccccccccccccccccccccccccccccc

  EPS=1.0E-10

  DEL=0.001

  JJ=0

  KOLSTE=30000

  H=0.001

  CALL START(XN, YN, VN, ZN)

C START VALUE SSSSSSSSS START VALUE

  YN7=YN+0.000001

          YN77=YN+0.000001

  YK=YN

  VK=VN

  ZK=ZN

  open(9,form='FORMATTED',file='RESU1.dat')

  open(8,form='FORMATTED',file='UGOL. dat')

  open(10,form='FORMATTED',file='PORTRET. dat')

  open(11,form='FORMATTED',file='XXXX. dat')

  open(6,form='FORMATTED',file='CKOROCT. dat')

  open(7,form='FORMATTED',file='rel1.for')

  WRITE (7,7) YK, VK, ZK 

  N=0

  10 N=N+1

  XNP1=XN+H

  K=0

  20 K=K+1

  DGDY=(G(XNP1,YK+DEL, VK, ZK)-G(XNP1,YK, VK, ZK))/DEL

  DGIDY=1.-H/2.*DGDY

  DGDV=(G(XNP1,YK, VK+DEL, ZK)-G(XNP1,YK, VK, ZK))/DEL

  DGIDV=-H/2.*DGDV

C

  DGDZ=(G(XNP1,YK, VK, ZK+DEL)-G(XNP1,YK, VK, ZK))/DEL

  DGIDZ=-H/2.*DGDZ

c

  DFDY=(F(XNP1,YK+DEL, VK, ZK)-F(XNP1,YK, VK, ZK))/DEL

  DFIDY=-H/2.*DFDY

c

  DFDV=(F(XNP1,YK, VK+DEL, ZK)-F(XNP1,YK, VK, ZK))/DEL

  DFIDV=1.-H/2.*DFDV

C

  DFDZ=(F(XNP1,YK, VK, ZK+DEL)-F(XNP1,YK, VK, ZK))/DEL

  DFIDZ=-H/2.*DFDZ

C

  DEDY=(E(XNP1,YK+DEL, VK, ZK)-E(XNP1,YK, VK, ZK))/DEL

  DEIDY=-H/2.*DEDY

c

  DEDV=(E(XNP1,YK, VK+DEL, ZK)-E(XNP1,YK, VK, ZK))/DEL

  DEIDV=-H/2.*DEDV

C

  DEDZ=(E(XNP1,YK, VK, ZK+DEL)-E(XNP1,YK, VK, ZK))/DEL

  DEIDZ=1.-H/2.*DEDZ

C

  DELTA=CA(DGIDY, DFIDY, DEIDY, DGIDV, DFIDV, DEIDV, DGIDZ, DFIDZ, DEIDZ)

C

  GIK=YK-H/2.*G(XNP1,YK, VK, ZK)-YN-H/2.*G(XN, YN, VN, ZN)

  FIK=VK-H/2.*F(XNP1,YK, VK, ZK)-VN-H/2.*F(XN, YN, VN, ZN)

  EIK=ZK-H/2.*E(XNP1,YK, VK, ZK)-ZN-H/2.*E(XN, YN, VN, ZN)

C

  DELTY=CA(-GIK,-FIK,-EIK, DGIDV, DFIDV, DEIDV, DGIDZ, DFIDZ, DEIDZ)

  DELTV=CA(DGIDY, DFIDY, DEIDY,-GIK,-FIK,-EIK, DGIDZ, DFIDZ, DEIDZ)

  DELTZ=CA(DGIDY, DFIDY, DEIDY, DGIDV, DFIDV, DEIDV,-GIK,-FIK,-EIK)

C

  YKP1=YK+DELTY/DELTA

  VKP1=VK+DELTV/DELTA

  ZKP1=ZK+DELTZ/DELTA

C

  ABSER=ABS(YKP1-YK)

  RELER=ABS(ABSER/YKP1)

  YK=YKP1

  VK=VKP1

  ZK=ZKP1

c  WRITE (9,7) DGIDY, DFIDY, DEDY, DGIDV, DFIDV, DEDV, DGIDZ, DFIDZ, DEDZ

  IF(RELER. GT. EPS. AND. K.LT.7)GOTO 20

ccccccccccccccccccccccccccccccccccccccc

       KS=SIGN(1.0,YN)

       KS1=SIGN(1.0,YKP1)

  IF(KS. NE. KS1.AND. KTURN. EQ.0)THEN

        SUM=0.0

               EE=(0.2/3*DD*DD+0.02*DD*DD)/2.0*VN*VN

                V2DV1=VKP1/VN

       KTURN=1

       KSW=0

               TIME=XNP1

       ENDIF

       KKS=SIGN(1.0,VN)

       KKS1=SIGN(1.0,VKP1)

        IF(KKS. NE. KKS1.AND. KSW. EQ.0)THEN

       DYDY=(YKP1-YN77)/YN77

        WRITE (7,7) XN, DYDY, YN77,YN7,YKP1,SUM, EE

                WRITE (*,7) XN, DYDY, YN77,YN7,YKP1,SUM, EE

       YN77=YN7

               YN7=YN

       KTURN=0

       KSW=1

        TIME=XNP1

       ENDIF

cccccccccccccccccccccccccccccccccccccccc

  XN=XNP1

  YN=YKP1

  VN=VKP1

  ZN=ZKP1

CC

ccccccccccccccccccccccccccccccccc

c  SUM=SUM+VM*VN**2*XX*H

  RR=DD-VV*(XNP1-TIME)

       IF((XNP1-TIME).LT. T0)THEN

  SUM=SUM+VM*VN**2*RR*H

        ENDIF

ccccccccccccccccccccccccccccccccc

  JJ=JJ+1 

  IF(JJ. EQ.5)THEN

c  WRITE (9,6) XN, YN, VN, ZN, K,N, ABSER

  WRITE (8,6) XN, YN

  WRITE (6,6) XN, VN

  WRITE (10,6) YN, VN

  WRITE (11,6) XN, XX

  JJ=0

  ENDIF

  IF(N. LE. KOLSTE)GOTO 10

  stop

  END

C 

          FUNCTION F(X, Y,V, Z)

  COMMON/CCC/KTURN, KSW, DD, ALP, ALM, TIME, BB, CC, CD3,AmM, XNP1,T0,XX, ALS

        T=XNP1-TIME

c                 CL=CD3+AmM*DD**2

               IF(T. GE. T0.AND. KTURN. EQ.1)THEN

        CL=CD3+AmM*BB**2

        XX=0.0

                KTURN=0

       AAA=0.0

         ENDIF

               

                       IF(T. LT. T0.AND. KTURN. EQ.1)THEN

       CL=ALP-ALS*T

        AAA=-1.0

                XX=(CL-CD3)/AmM

        ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccc

                       IF(T. GE. T0.AND. KSW. EQ.1)THEN

        CL=CD3+AmM*DD**2

        XX=0.0

                KSW=0

                CL=CD3+AmM*DD**2

                AAA=0.0

         ENDIF

               

                       IF(T. LT. T0.AND. KSW. EQ.1)THEN

       CL=ALM+ALS*T

        AAA=0.0

                CL=CD3+AmM*DD**2

        ENDIF

c                 CL=CD3+AmM*DD**2

       XX=CL

  F=-SIN(Y)/(CL*CC)-V*ALS/CL*AAA

  RETURN

  END

C

  FUNCTION G(X, Y,V, Z)

  G=V 

  RETURN

  END

C

  FUNCTION E(X, Y,V, Z)

  E=Y+V-Z-X+1. 

  RETURN

  END

c

  SUBROUTINE START(X, Y,V, Z)

  X=0.

c  Y=3.1414926/18.0  *2.5

  Y=3.1414926/18.0

  V=0.

  Z=2.

  RETURN

  END

C

  FUNCTION CA(A1,A2,A3,B1,B2,B3,C1,C2,C3)

  RAB1=A1*B2*C3+A2*B3*C1+A3*B1*C2

  RAB2=A3*B2*C1+A2*B1*C3+A1*B3*C2

  CA=RAB1-RAB2

  RETURN

  END