18. МАТЕМАТИКА

Язык программирования ACL является языком, с одной стороны, интерпретируемым, а с другой -- в известной степени специальным, то есть направленным на определенный класс задач, которые его автору приходится решать относительно часто. И хотя на ACL можно написать программу на любые вычисления, тем не менее есть резон какие-то вычисления общего порядка выполнять запуская готовые программы с помощью команд. В основном, это касается вычислений в длинных циклах, но часто речь может идти и просто о вычислениях специального вида. Использование команд не только ускоряет работу программы, но и делает код более читаемым. Все готовые вычисления реализуются в рамках одной команды называемой #mathematics, а короче #ma. Эта команда может иметь много операций и число ее операций растет от версии к версии. В этом разделе рассмотрены все операции, разработанные для последней версии программы. Операции можно условно разделить на группы. Первая группа операций -- это векторные операции. Под вектором будем понимать часть вещественного массива, начиная с индекса, задаваемого параметром [b] и с числом элементов [le]. Когда в операции используются и другие векторы, то они описываются другими параметрами. Название операции обычно состоит из первых букв английских слов описания операции. Полная расшифровка дается в описании каждой операции. Итак рассмотрим векторные операции.

.

      ВЕКТОРНЫЕ ОПЕРАЦИИ

.

 #ma [op=vic; b=; le=;]

(vector initialize constant). Эта операция выполняет инициализацию вектора константой. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Постоянное значение должно быть задано в переменной [C]. В результате операции все элементы выбранной части массива r() будут равны [C].

.

 #ma [op=via; b=; le=;]

(vector initialize arithmetic). Эта операция выполняет инициализацию вектора арифметической прогрессией. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Первое значение задается в переменной [C]. Шаг к следующему значению задается в переменной [D]. В результате операции все элементы выбранной части массива r() будут равны C + D*i, где i есть разность индексов текущего и первого элементов, то есть он равен 0 для первого элемента.

.

 #ma [op=vig; b=; le=;]

(vector initialize Gauss). Эта операция выполняет инициализацию вектора функцией Гаусса. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Первое значение аргумента задается в переменной [A]. Шаг к следующему значению аргумента задается в переменной [B]. Положение центра кривой задается в переменной [C]. Полуширина, то есть ширина на половине высоты задается в переменной [D]. В результате операции все элементы выбранной части массива r() будут равны значениям функции exp(-a1*(A+B*i-C)^2), где i есть разность индексов текущего и первого элементов, то есть он равен 0 для первого элемента, a1 - коэффициент, вычисляемый из полуширины кривой D.

.

 #ma [op=vaс; b=; le=;]

(vector add constant). Эта операция выполняет прибавление к вектору константы. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Постоянное значение должно быть задано в переменной [C]. В результате операции ко всем элементам выбранной части массива r() добавляется константа [C].

.

 #ma [op=vmс; b=; le=;]

(vector multiply constant). Эта операция выполняет умножение вектора на константу. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Постоянное значение должно быть задано в переменной [C]. В результате операции все элементы выбранной части массива r() умножаются на константу [C].

.

 #ma [op=vex; b=; le=;]

(vector exponent). Эта операция вычисляет экспоненту от значений вектора поэлементно. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. В результате операции все элементы выбранной части массива r() получают новые значения.

.

 #ma [op=vlo; b=; le=;]

(vector logarithm). Эта операция вычисляет натуральный логарифм от значений вектора поэлементно. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. В результате операции все элементы выбранной части массива r() получают новые значения.

.

 #ma [op=vsi; b=; le=;]

(vector sine). Эта операция вычисляет синус от значений вектора поэлементно. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. В результате операции все элементы выбранной части массива r() получают новые значения.

.

 #ma [op=vco; b=; le=;]

(vector cosine). Эта операция вычисляет косинус от значений вектора поэлементно. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. В результате операции все элементы выбранной части массива r() получают новые значения.

.

 #ma [op=vpo; b=; le=;]

(vector power). Эта операция вычисляет степень от значений вектора поэлементно. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Показатель степени должен находиться в переменной C. Проверка на неправильные операции не проводится. Пользователь должен заботиться об этом сам. В результате операции все элементы выбранной части массива r() получают новые значения, то есть x превращается в x^C.

.

 #ma [op=vin; b=; le=;]

(vector integer). Эта операция вычисляет целую часть снизу от значений вектора поэлементно. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. В результате операции все элементы выбранной части массива r() получают новые значения. Эта операция необходима потому что пересылка реального массива в целый дает ближайшее целое. А иногда нужно знать именно целую часть снизу.

.

 #ma [op=vba; b=; le=;]

(vector byte ASCII). Эта специальная операция выполняет преобразование вектора из байтов в ASCII коды. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Предполагается, что в нем записаны байты, то есть целые числа в диапазоне от -128 до 127. В результате преобразования отрицательные числа заменяются на положительные прибавлением 256.

.

 #ma [op=vva; b=; le=; trx=;]

(vector vector addition). Эта операция выполняет сложение вектора с вектором поэлементно. Первый вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Второй вектор такой же длины смещен относительно первого на разность индексов [trx]. Результат сложения записывается вместо первого вектора. Параметр [trx] может быть отрицателен. Для правильной работы векторы не должны перекрываться. Для вычитания векторов достаточно второй вектор предварительно умножить на -1.

.

 #ma [op=vvm; b=; le=; trx=;]

(vector vector multiply). Эта операция выполняет умножение вектора на вектор поэлементно. Первый вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Второй вектор такой же длины смещен относительно первого на разность индексов [trx]. Результат умножения записывается вместо первого вектора. Параметр [trx] может быть отрицателен. Для правильной работы векторы не должны перекрываться, но возможно когда они совпадают, то есть trx=0. В этом случае вычисляется квадраты всех элементов.

.

 #ma [op=vvs; b=; le=; trx=;]

(vector vector scalar). Эта операция выполняет скалярное умножение вектора на вектор. Первый вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Второй вектор такой же длины смещен относительно первого на разность индексов [trx]. Результат умножения записывается в переменную [C].

.

 #ma [op=var; b=; le=;]

(vector area). Эта операция ограничивает область значений вектора. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Если значение элементов вектора меньше, чем значение переменной V0, то оно заменяется на V0. Соответственно, если больше V1, то заменяется на V1.

.

 #ma [op=vor; b=; le=; mo=;]

(vector order). Эта операция упорядочивает элементы вектора по возрастанию, если [mo] равно 1, или убыванию, если [mo] равно 0. Вектор находится в реальном массиве r() и определяется параметрами [b] и [le]. Одновременно в целом массиве i() начиная с первого элемента указываются индексы элементов, какие они имели до упорядочивания. Таким образом, если необходимо упорядочить аналогичным образом другие векторы, то это уже можно сделать используя массив i(). При этом первым надо поставить тот индекс, который находится в i(1) и так далее.

.

 #ma [op=cwm; b=; le=; mo=;]

(curve width maximum). Эта операция выполняет расчет параметров максимума (фокуса) кривой, заданной некоторой системой значений реального массива r(). Параметры есть [b] как индекс первого элемента массива r(), и [le] как число рассматриваемых значений. В результате операции будут вычислены минимальное и максимальное значения, центр HM-области, и ширина HM-области. Эти значения будут записаны в переменные A,B,C,D. HM-область (half-maximum) -- это область, внутри которой значения элементов массива больше чем B/2 если параметр [mo] равен 0 или B-(B-A)/2, если [mo] равен 1. Максимальное значение вычисляется интерполяцией кривой по ближайшим точкам. Оно всегда больше любого из реальных значений этого массива. В этом отличие этой операции от операции (op=fno), где просто выбирается максимальное значение из значений элементов массива.

.

 #ma [op=fsd; b=; le=; trx=; n=;]

(Function from Second Derivative). Эта операция вычисляет функцию по двум точкам и значениям ее второй производной. Сама функция находится в массиве r() начиная с индекса [b] и длиной [le=N;]. Две точки функции должны быть заданы начиная с индекса [n], индекс отсчитывается от начала массива функции и имеет первое значение 1, то есть индексу массива r(), заданному в [b], соответствует n=1. Остальные точки вычисляются по формуле f[k+1] = F[k] + 2*f[k] -- f[k-1], k=n+1,...,N-1 и f[k-1] = F[k] + 2*f[k] -- f[k+1], k=n,...,2. Вектор функции F[k] должен быть задан в массиве r(), начиная с индекса [b]+[trx]. Если n=1, то функция вычисляется по первым двум точкам слева направо. Если n > 1, то сначала вычисляются все правые точки, а затем левые. Этот случай необходим, если на краях интервала вторая производная определена неправильно. Тогда движение от края приводит к неправильному результату во всех точках.

.

 #ma [op=fdf; b=; le=; trx=;]

(First Derivative of Function). Эта операция вычисляет первую производную от функции. Сама функция находится в массиве r() начиная с индекса [b]+[trx] и длиной [le=N;]. Первая производная вычисляется по формуле f[k] = (F[k+1] -- F[k-1])/2, k=2,...,N-1. Затем f[1]=f[2], f[N]=f[N-1]. Первая производная функции f[k] определяется в массиве r(), начиная с индекса [b]. Шаг массива точек считается равным 1, если это не так, то после вычисления на него надо делить весь массив.

.

 #ma [op=sdf; b=; le=; trx=;]

(Second Derivative of Function). Эта операция вычисляет вторую производную от функции. Сама функция находится в массиве r() начиная с индекса [b]+[trx] и длиной [le=N;]. Вторая производная вычисляется по формуле f[k] = F[k+1] + F[k-1]-2*F[k], k=2,...,N-1. Затем f[1]=f[2], f[N]=f[N-1]. Вторая производная функции f[k] определяется в массиве r(), начиная с индекса [b]. Шаг массива точек считается равным 1, если это не так, то после вычисления надо делить весь массив на квадрат шага.

      МАТРИЧНЫЕ ОПЕРАЦИИ

Вторая группа операций -- это матричные операции. Под матрицей понимается часть вещественного массива r() с начальным индексом, задаваемым параметров [b], но размерность матрицы задается параметрами [nx] (число колонок) и [ny] (число строк).

.

 #ma [op=mtr; b=; nx=; ny=;]

(matrix transpose). Эта операция выполняет транспонирование матрицы. Матрица находится в реальном массиве r() с первого номера, задаваемого параметром [b]. Число строк задается [ny] а число столбцов [nx]. В результате операции в той же части реального массива будет записана транспонированная матрица у которой колонки станут строками, а строки колонками. Эта операция очень удобна при чтении чисел из файла в свободном формате (из колонок) в том случае, когда необходимо выделить отдельные функции из многих функций. После транспонирования значения каждой функции записаны слитно.

.

 #ma [op=mii; b=; nx=; ny=; mo=;]

(matrix index inversion). Эта операция выполняет инверсию индексов матрицы. В результате операции последний элемент каждой строки становится первым и последняя строка становится первой. Как обычно, матрица описывается параметрами [nx], [ny] и [b]. Кроме того, параметр [mo] определяет модификации. Если [mo=0;], оба индекса инвертируются, если [mo=1;], то только первый индекс инвертируется, то есть строки, если [mo=2;], то только второй индекс инвертируется, то есть колонки. Результирующая матрица записывается на то же место. Заметим, что предполагается обычный в математике порядок записи матрицы строками.

.

 #ma [op=mss; b=; nx=; ny=;]

(matrix sum sum). Эта операция выполняет суммирование строк и колонок матрицы, которая имеет форму m[nx,ny]. В результате вычисляется вектор v[nx] как сумма по всем k для m[nx,k] который записывается сразу за матрицей, затем вычисляется вектор v[ny] как сумма по всем k для m[k,ny] который записывается после предыдущего вектора. А следующий элемент реального массива равен полной сумме всех элементов матрицы. Это можно себе представить как интегрирование функции f(x,y) сначала по y с результатом A(x), затем по x с результатом B(y) и затем по x и y одновременно. В важном частном случае, если [ny] = 1 и матрица превращается в вектор, операция работает иначе. В этом случае вычисляется только сумма по всем элементам вектора и результат записывается за вектором. Это эквивалентно интегралу от функции f(x). Так это работает при значении параметра [emp=0;]. Если [emp=1;] то выполняется суммирование только по второму индексу при условии, что ny > 1. А при [emp=2;] выполняется суммирование только по первому индексу. Каждый раз результат суммирования записывается сразу за матрицей а за ним полная сумма. Новые режимы позволяют упростить программирование и ускорить счет в тех случаях, когда двойное суммирование не является необходимым.

.

 #ma [op=mmi; n=; b=; nx=; ny=; sav=; wid=; hei=; file=; form=; em=;]

(matrix matrix interpolation). Эта операция выполняет интерполяцию матрицы в матрицу. Это довольно сложная операция, потому что размер матрицы может быть произвольно большим. Начальная матрица имеет размеры [nx] и [ny], в то время как новая матрица имеет размеры [wid] и [hei]. Предполагается, что исходная матрица описывает прямоугольную область с xi1, yi1, xi2, yi2, в то время как интерполированная матрица описывает прямоугольную область xo1, yo1, xo2, yo2. Эти значения должны быть заданы в реальном массиве r() начиная с индекса задаваемого параметром [n] в следующем порядке xi1,xi2,yi1,yi2,xo1,xo2,yo1,yo2. Обычно исходная матрица берется из файла с именем [file]. Однако в случае [file=here;] матрица должна быть размещена в реальном массиве r() стартуя с индекса [b]. Обычно результирующая матрица записывается в файл с именем [form] . Однако в случае [form=here;] матрица записывается в реальный массив r() стартуя с индекса [sav]. Предполагается что шаг сетки точек задается по формуле dx=(x2-x1)/(nx-1). Используется линейная интерполяция. Если размер области интеполированной матрицы превышают размеры области исходной матрицы, то значения в точках за пределами исходной матрицы задаются как средние значения исходной матрицы если параметр [em] не равен 1. А если равен, то они задаются равными значению переменной C.

.

 #ma [op=mmr; b=; le=;]

(matrix matrix rotation). Эта операция выполняет вращение квадратной матрицы вокруг ее центра. Предполагается, что квадратная матрица записана в реальном массиве r() начиная с индекса [b]. Размерность матрицы равна [le]*[le]. Из этой матрицы получается другая матрица, такой же размерности, элементы которой вычисляются из элементов исходной матрицы по закону
xo=xn*cos(A)-yn*sin(A), yo=xn*sin(A)+yn*cos(A).
Здесь угол A (в радианах) равен значению переменной [A], xn,yn - координаты новой матрицы, то есть индексы элемента, отсчитываемые от середины. Если [le] четное, то эти координаты полуцелые. Вычисленные координаты xo,yo старой матрицы используются для определения индексов 4-х ближайших элементов, из которых производится интерполяция. Если индексы выходят за пределы их изменения, то соответствующие элементы полагаются равными нулю. Операция эквивалентна вращению черно-белой картинки, которую задает матрица. Новая матрица записывается в том же массиве сразу за исходной, то есть начиная с элемента [b]+[le]*[le].

.

 #ma [op=mcg; b=; sav=; nx=; ny=; n=i; file=; form=;]

(matrix convolution Gauss). Эта операция выполняет свертку (или развертку) матрицы с гауссианом. Предполагается, что матрица используется в вычислениях типа FFT (быстрого преобразования Фурье) и она имеет размерности [nx] и [ny] как целые степени 2, то есть 64, 128, 256, 512, 1024. Однако [ny] может быть равно 1 (см. далее). Светка (или развертка) вычисляется методом двойного FFT и результатом является матрица таких же размеров. Параметры гауссиана (sigmaX) и (sigmaY) измеряются в единицах равных размеру одного пиксела изображения. Гауссиан задается по каждому измерению своим преобразованием Фурье
G(q) = exp( -(sigma*q)^2/2 ).
Свертка выполняется при условии, если параметр (epsilon) = 0. В этом случае фурье-изображение матрицы умножается на G(q). В противном случае если (epsilon) > 0, выполняется развертка при разумной регуляризации, а именно, Фурье изображение матрицы умножается на множитель
1 / sqrt( G(q)^2 + epsilon^2 ).
Значения (sigmaX), (sigmaY), (epsilon) должны быть заданы в элементах реального массива r(i), r(i+1) и r(i+2),, где (i) определяется параметром [n]. Замечу, что если (sigmaX) = 0 или (sigmaY) = 0, то по соответствующему направлению ничего не делается. Исходная матрица читается из файла с именем [file]. Если [file=here;], матрица читается из реального массива начиная с индекса [b]. Вычисленная матрица записывается в файл, имя которого задается параметром [form]. Если [form=here;] матрица будет записана в реальный массив r() стартуя с индекса [sav]. Значение [ny=1;] также возможно и означает одномерный случай, когда матрица совпадает с вектором размера [nx]. Исходный файл должен быть записан в режиме float, то есть 4 байта на число. Результат записывается в файл в таком же режиме.

.

 # C=p1; D=p2; 
 #ma [op=mip; b=; nx=; ny=;]

(matrix initialize periodicity). Эта операция выполняет преобразование матрицы, которое условно можно назвать размножением периода. Предполагается, что матрица записана в массив r(), начиная с индекса [b] и имеет размеры [nx] и [ny] по горизонтали и вертикали. Необходимо выделить центральную часть этой матрицы внутри прямоугольника с размерами p1 (по горизонтали) и p2 (по вертикали) и переопределить все элементы матрицы вне этого прямоугольника из условия периодичности. То есть сдвигать координаты точки на целое число периодов, до тех пор, пока они не попадут в центральный прямоугольник и потом определить значение точки внутри прямоугольника интерполяцией. Сначала определяются области по вертикали, потом по горизонтали. Периоды в общем случае могут не равняться целому числу шагов матрицы, они задаются в переменных C и D в единицах, в которых один элемент матрицы имеет размер 1. Применение такой процедуре можно найти в разных задачах. Например, пусть вам надо задать периодическую картинку в большой матрице. Тогда достаточно просто исходно задать нулевую матрицу, определить один период и использовать процедуру. Или у вас уже задана периодическая матрица, но в результате ее преобразований периодичность по краям пропала. Эта процедура позволяет исправить ситуацию. В графике аналогичная процедура называется заполнением области текстурой. Задается произвольная область, задается картинка и картинка копируется в области периодически начиная с левого верхнего угла. Мне было необходимо копировать именно из центра области. Чтобы после операции матрица не имела скачков необходимо, чтобы в центральной части значения функции на границах периода совпадали. Одномерный случай тоже нормально работает при условии, что задаются значения ny=1, D=1.

.

 #ma [op=fno; b=; nx=; ny=; mo=;]

(function normalization). Эта операция выполняет специальную нормировку матрицы. Матрица определяется параметрами [b], [nx], [ny], как и выше. Однако в этом случае каждая строка размером [nx] трактуется как вектор v(j). Определяются минимальное v_min и максимальное v_max значения этого вектора и [ny] значений v_min и [ny] значений v_max записываются сразу за матрицей, начиная с индекса [b]+[nx]*[ny]. Затем если [mo=0;], то каждый вектор нормируется согласно формуле [v(j) - v_min]/[v_max - v_min]. Нормированные значения могут быть использованы для графики. Если [mo] не равно нулю, то матрица остается той же самой и операция используется только для определения минимальных и максимальных значений.

.

 #ma [op=ffi; b=; nx=; ny=; n=; wid=; c=;]

(function function interpolation). Эта операция выполняет интерполяцию строк матрицы, описывающих функции на неравномерной сетке в строки матрицы, описывающие функции на равномерной сетке и в произвольном интервале с любым числом точек, включая одну точку. Интерполяция линейная в пределах интервала содержащего новую точку. Если новая точка находится за пределами области определения функций, то берется соответствующее граничное значение. Предполагается, что исходная матрица записана в реальном массиве r() начиная с индекса [b] и имеет размеры [nx] и [ny], в то время как новая матрица имеет размеры [wid] и [ny] и записывается в r(), начиная с индекса [n]. Строка аргументов исходной матрицы записана в r() непосредственно перед ней, а аргументы новой матрицы определяются по начальному и конечному значениям, которые должны быть заданы в r() начиная с индекса, задаваемого параметром [c]. Эта операция удобна при конвертировании внешних таблиц в таблицы, используемые в расчете.

.

 #ma [op=fis; b=; nx=; ny=; siz=;]

(function interpolate spline). Эта операция выполняет интерполяцию строк матрицы кубическим сплайном. Каждая строка исходной матрицы представляет собой значения функции, определенной на системе значений аргумента с переменным шагом. Матрица имеет размеры [nx] и [ny] и записана в реальном массиве r() начиная с индекса [b]. Аргументы должны предшествовать матрице и система точек начинается с индекса [b]-[nx]. Результатом интерполяции является новая матрица, которая имеет размеры [siz] and [ny]. Она определяет ту же самую систему функций, но на новой системе точек с постоянным шагом начиная с первой старой точки до последней старой точки. Новая матрица записывается сразу за старой матрицей, так что для нее первый элемент массива r() начинается с [b]+[nx]*[ny]. Другими словами, эта операция меняет представление системы функций с переменного шага на постоянный шаг с использованием кубического сплайна.

.

 #ma [op=mds; b=; nx=; ny=; xp=; yp=; sa=; n=; dir=;]

(matrix diagonal selection). Эта операция выполняет специальный расчет, который состоит в следующем: Выделяется одномерный массив с началом в [sa] и размером [n]. Первоначально этот массив обнуляется. Затем в середину этого массива записывается значение элемента матрицы с индексами [xp] и [yp]. Индексы начинаются от единицы. Сама матрица должна быть заранее записана с началом в [b] и размерами [nx] и [ny]. Затем все точки одномерного массива заполняются элементами матрицы, находящимися на диагонали, проходящей через указанную точку до тех пор, пока диагональ не выйдет за размеры матрицы. То есть заполнение идет вправо от центральной точки до выхода за границы массива и затем влево таким же способом. При этом одномерный массив заполняется не полностью, и часть его элементов остаются нулями. При dir=1 диагональ идет вправо и наверх, то есть (1,1), а при dir=-1 она идет влево и наверх, то есть (-1,1).Операция полезна когда матрица описывает локальную функцию и ее необходимо проинтегрировать по диагонали с весом.

.

 #ma [op=smh; file=; le=; twi=; the=; ord=; b=;]

(search maximum hexadecimal). Эта операция выполняет поиск области максимума (фокуса) внутри 16-битного изображения, записанного детектором. Это достаточно специальная операция, но она может быть полезной в некоторых случаях. Параметры означают следующее: [file] -- имя файла, [le] -- размер заголовка в байтах, [twi],[the] -- размеры изображения в пикселах, [ord] -- порядок байтов в двухбайтовом представлении чисел, [ord=0;] для Low Byte Second и [ord=1;] для Low Byte First. Результат будет записан в реальный массив r(), начиная с индекса [b]. Результат представляет собой 5 чисел: x и y координаты фокуса, значение максимума в фокусе, а также x и y компоненты FWHM (полная ширина на половине высоты) для пика фокуса. Все координаты измеряются в пикселах.

.

 #ma [op=rfh; file=; le=; twi=; the=; ord=; b=; wid=; hei=; xs=; ys=;]

(read fragment hexadecimal). Эта операция выполняет чтение фрагмента внутри 16-битного изображения, записанного детектором. Это достаточно специальная операция, но она может быть полезной в некоторых случаях. Параметры означают следующее: [file] -- имя файла, [le] -- размер заголовка в байтах, [twi],[the] -- размеры изображения в пикселах, [ord] -- порядок байтов в двухбайтовом представлении чисел, [ord=0;] для Low Byte Second и [ord=1;] для Low Byte First. Результат будет записан в реальный массив r(), начиная с индекса [b]. Результат представляет собой прямоугольную матрицу чисел шириной [wid] и высотой [hei]. Если трактовать исходную картинку как записанную снизу вверх, то [xs] и [ys] задают координаты левой нижней точки фрагмента, начальные значения равны (1,1). Если ширина или высота фрагмента равны единице, то результатом будет часть столбца или строки исходной матрицы.

      ОПЕРАЦИИ С КОМПЛЕКСНЫМИ ЧИСЛАМИ

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

.

 #ma [op=mcc; b=; nx=;]

(multiply complex complex). Эта операция умножает комплексный массив на комплексный массив. Первый комплексный массив находится в реальном массиве r(), начиная с элемента с индексом [b]. Он имеет структуру r,i,r,i,r,i и размер [nx] комплексных значений. Второй комплексный массив расположен сразу после первого. Результат умножения будет помещен на место первого массива, так что первый массив будет потерян. Это устаревшая команда, которая сохранена для работы старых ACL программ. Вместо нее в версиях после 6.4.13 можно использовать команду

.

 #ma [op=ccm; b=; le=; trx=;]

(complex complex multiplication). Эта операция умножает комплексный массив на комплексный массив. Первый комплексный массив находится в реальном массиве r(), начиная с элемента с индексом [b]. Он имеет структуру r,i,r,i,r,i и размер [le] комплексных значений. Второй комплексный массив расположен со сдвигом индекса на [trx]. Результат умножения будет помещен на место первого массива, так что первый массив будет потерян. Обратите внимание: при переписывании старой команды на новую надо размер массива указывать в [le] вместо [nx].

.

 #ma [op=dcc; b=; nx=;]

(division complex complex). Эта операция делит комплексный массив на комплексный массив. Первый комплексный массив находится реальном массиве r() начиная с элемента с индексом [b]. Он имеет структуру r,i,r,i,r,i и размер [nx] комплексных значений. Второй комплексный массив расположен сразу после первого. Результат деления будет помещен на место первого массива, так что первый массив будет потерян в памяти. Чтобы избежать деления на ноль операция делается следующим образом c1*conjg(c2)/(abs(c2)^2+E^2). Поэтому должна быть определена переменная [E], которая должна иметь достаточно малое значение. Это устаревшая команда, которая сохранена для работы старых ACL программ. Вместо нее в версии 6.4.13(140) можно использовать команду

.

 #ma [op=ccd; b=; le=; trx=;]

(complex complex division). Эта операция делит комплексный массив на комплексный массив. Первый комплексный массив находится реальном массиве r() начиная с элемента с индексом [b]. Он имеет структуру r,i,r,i,r,i и размер [le] комплексных значений. Второй комплексный массив расположен со сдвигом индекса на [trx]. Результат деления будет помещен на место первого массива, так что первый массив будет потерян в памяти. Чтобы избежать деления на ноль операция делается следующим образом c1*conjg(c2)/(abs(c2)^2+E^2). Поэтому должна быть определена переменная [E], которая должна иметь достаточно малое значение. Обратите внимание: при переписывании старой команды на новую надо размер массива указывать в [le] вместо [nx].

.

 #ma [op=cen; b=; n=; trx=;]

(conversion exponent normal). Эта операция преобразует комплексный массив чисел из экспоненциальной в нормальную форму. Размер массива определяется параметром [num]. В экспоненциальной форме каждый элемент массива задается модулем и фазой, причем в реальном массиве r() сначала записывается [n] значений модуля, начиная с индекса, определяемого параметром [beg], а сразу вслед за ним [n] значений фазы для тех же комплексных чисел. Эти данные должны быть заданы до начала операции. В результате операции в реальный массив r() записывается [n] комплексных чисел в режиме r,i,r,i,... причем первый элемент записывается в индекс [b]+[trx]. Очевидно старая и новая записи не должны перекрываться, поэтому [trx] по модулю должно быть больше 2*[n]. Однако возможная ощибка не отслеживается, будьте внимательны. Данная операция удобна, если вам необходимо задавать комплексный массив через зависимость модуля или фазы от какого-либо параметра. Так как в расчетах фурье-преобразования используется нормальная форма, то такая операция необходима.

.

 #ma [op=cne; b=; n=; nx=; ny=; trx=;]

(conversion normal exponent). Эта операция является обратной к предыдущей и она преобразует комплексный массив чисел из нормальной в экспоненциальную форму. Размер массива определяется параметром [n]. В экспоненциальной форме каждый элемент массива задается модулем и фазой, причем в реальном массиве r() сначала записывается [n] значений модуля, начиная с индекса, определяемого параметром [b]+[trx], а сразу вслед за ним [n] значений фазы для тех же комплексных чисел. Эти данные являются результатом операции. А до начала операции в реальный массив r() записывается [n] комплексных чисел в режиме r,i,r,i,... причем первый элемент записывается в индекс [b]. Очевидно старая и новая записи не должны перекрываться, поэтому [trx] по модулю должно быть больше 2*[n]. Однако возможная ощибка не отслеживается, будьте внимательны. Данная операция удобна, если вам необходимо проанализировать графически комплексный массив через зависимость модуля или фазы от какого-либо параметра. Так как в расчетах фурье-преобразования используется нормальная форма, то такая операция необходима. Для удобства графического представления, а также возможно и для других целей массив реальных чисел, описывающих фазу выпрямляется, то есть ликвидируются возможные скачки на 2 пи. Эти скачки неизбежны, так как фаза определяется через arctan(i/r) с учетом конкретных знаков r и i. Алгоритм выпрямления фазы зависит от параметра [nx]. Значение этого параметра должно находиться между 1 и [n] а его смысл в том, что он задает стартовую точку для выпрямления фазы как в сторону возрастания индекса, так и в сторону убывания. Эта стартовая точка должна находиться в области предполагаемого наименьшего изменения фазы. Например, для симметричной параболической фазы она находится в середине области, то есть [nx]=[n]/2. Это значение для [nx] автоматически устанавливается без указания на ошибку, если исходное значение [nx] находится вне указанного интервала. При конвертировании двумерного массива механизм выпрямления фазы нужно отключать. Он отключается, если [ny]=10101. Случайно такое значение поставить маловероятно, поэтому можно не заботиться если не надо отключать выпрямление. Но если надо, то его надо ставить принудительно.

.

 #ma [op=fft; b=; nx=; em=;]

(fast Fourier transformation). Эта операция означает Фурье преобразование в бесконечных пределах для функции, которая равна нулю на бесконечности. Расчет проводится с использованием алгоритма FFT (быстрого преобразования Фурье) для суммы
fq[n] = sum{m = 0,...,N-1} exp( sign(k)*2*pi*i*n*m/N )*fx[m]
где N = 2^abs(k) и k -- целое, n = 0,...,N-1
Число k = log(N)/log(2) задается в параметре [nx]. При этом k может быть отрицательным, но его абсолютное значение не может быть больше 20. При k < 0 вычисляется переход из x пространства в q пространство. При этом функция в x-пространстве должна быть задана на системе точек x[n]=D*(n+(1-N)/2), где D - шаг сетки точек - задается в программе в указанной переменной [D]. При этом функция определяется в симметричных пределах на области X = D*N, а для преобразования интеграла к сумме используются некоторые дополнительные умножения на фазовые множители. Преобразования не зависят от значения D и оно необходимо только для умножения на сумму и получения ответа в физической размерной системе координат. Ответ получается на системе точек q[n]=d*(n+(1-N)/2), где d = 2п/X, а полная область в q пространстве равна Q = 2п/D. Если k > 0, то все наоборот. Функция должна быть задана на системе точек q[n]=d*(n+(1-N)/2), где d = 2п/X, ответ получается на системе точек x[n]=D*(n+(1-N)/2), где D = X/N, а сумму надо умножать на множитель d/2п = 1/X. Этот множитель также следует задать в переменной [D]. То есть переменная [D] есть параметр операции. Комплексный массив fx[m] должен быть записан в реальный массив r(), начиная с элемента с индексом [b] и упорядоченный как (r,i,r,i, и так далее). Вычисленное фурье-изображение будет помещено к той же части реального массива r().

Процедура имеет несколько вариантов, которые различаются значениями параметра [em]. Если [em=0], то вычисляется интеграл по старому коду программы с использованием переменных комплексного типа. При расчете интеграла массив функции до после вычисления суммы умножается на массив экспоненциальных множителей и результат получается правильный. Если [em=1], то вычисляется только указанная выше сумма. При этом массив результата получается сдвинутым на полпериода и искаженным таким образом, что каждое следующее значение меняет знак. Но обратное преобразование точно дает исходную функцию. Расчет суммы и является процедурой FFT и, например, в Питоне только таким расчетом и ограничиваются. Это вариант сделан как раз для сравнения с другими языками. Если [em=2], то снова вычисляется интеграл, но без переменных комплексного типа. Расчет дает такой же ответ, но чуть быстрее по времени. Этот вариант был сделан в 2020 году. Если [em=3], то снова вычисляется сумма, но без переменных комплексного типа. При этом расчет выполняется в два раза быстрее, чем при [em=1. Это вариант был сделан одновременно с предыдущим вариантом.

.

 #ma [op=com; em=0;]  K A B C   K A B C   K A B C   K A B C . . .

(complex). Это операция программирования арифметики комплексных чисел. Она задает порядок операций над комплексными числами методом записи команд как было в первых компьютерах. Каждая команда содержит 4 целых числа K A B C, причем первое число задает код операции, а остальные три - указания на размещение аргументов в реальном массиве r(), то есть на три комплексных числа, причем над первыми двумя производится операция, а результат записывается в третье число. Другими словами, комплексные числа записаны в реальном массиве r(), а целые числа A, B, C - это индексы для реальной части рассматриваемых комплексных чисел, мнимые части следуют за ними. Реализованы следующие коды операций: 1 - сложение, 2 - вычитание, 3 - умножение, 4 - деление, 5 - exp(A), 6 - cos(A), 7 - sin(A), 8 - tan(A), 9 - pow(A,B), 10 - log(A), 11 - J0(A), 12 - J1(A). Для функций второй аргумент не используется. Исключением является функция pow(), для которой используется реальная часть второго аргумента, то есть комплексное число возводится только в реальную степень. J0 и J1 -- это функции Бесселя нулевого и первого порядков.

Вот как это работает. Пусть вам необходимо вычислить выражение a=a+b*exp(c)/(d+e); где a,b,c,d,e - комплексные числа. Вы вводите еще два числа f и g и все числа размещаете в реальном массиве так что a=(r(1),r(2)), b=(r(3),r(4)), . . ., g=(r(13),r(14)). Тогда выражение можно вычислить следующей системой команд

 #ma [op=com; em=0;] 5 5 0 11   3 11 3 11   1 7 9 13   4 11 13 11   1 1 11 1

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

 #ma [op=com; em=1; b=;]  + a b c   * a c c   S a a c   0 a a c . . .

Старый вариант стал работать при параметре [em=0;], а новый при [em=1;], и в новом варианте надо еще задать параметр [b] как первый номер массива комплексных чисел в реальном массиве. А сами комплексные числа записываются маленькими латинскими буквами, всего 26 от a до z. Соответствие с реальным массивом r() такое. Буква a имеет индекс, указываемый параметром [b]. Для буквы b индекс на 2 больше. Для с еще на 2 больше и так далее. Какие-то буквы надо задать до выполнения команды, какие-то будут вычислены после ее определения. Значения потом можно взять из реального массива. Что касается операций, то они вместо номеров задаются следующими символами:
+   -   *   /   E   C   S   T   P   L   0   1
Такие символы более наглядно указывают на операции и функции, а смысл у них тот же самый, какой указан выше для номеров. Например, E = exp(), L = log(), последние два символа указывают на функции Бесселя.

К сожалению, комплексная математика в ACL, как и в Java была развита слабо и позднее всего остального кода. Раньше я такие вычисления делал непосредственно на Java в виде готовых процедур. Но с течением времени ситуация упрощается. После мая 2021 года появилась команда скалярного умножения комплексных матриц. Она реализована в этой же операции, которую мы сейчас обсуждаем, но при параметре [em=2;]. Выглядит она очень просто.

 #ma [op=com; em=2; b=;]

Но это потому, что ее аргументы надо записывать в целый массив. То есть параметр [b=;] указывает на первый элемент целого массива i(), в котором надо записать 6 целых чисел: i1 i2 i3 n1 n2 n3. Они означают следующее. i1 -- это начальный индекс реального массива r(), где записана первая матрица. i2 -- то же самое для второй матрицы. i3 -- адрес реального массива, где будет записан ответ. Матрицы и ответ не должны пересекаться. Далее n1 -- это число строк первой матрицы, n2 -- число столбцов первой матрицы, вторая матрица имеет точно такое же число строк, наконец, n3 -- число столбцов второй матрицы. В частном случае эта операция позволяет умножать комплексный вектор на комплексное число, если размерности равны n 1 1. Или можно сделать матрицу из двух векторов при n 1 m.

      СПЕЦИАЛЬНЫЕ ВЫЧИСЛЕНИЯ

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

Специальные вычисления не являются элементами языка программирования в буквальном смысле этого слова. Все такие вычисления можно сделать с помощью других команд ACL. Но часто это не удобно по разным причинам. Одной из причин может быть массовость таких вычислений. Другой -- необходимость выполнить вычисления быстро. В связи с этим данная группа операций тоже постепенно увеличивается в числе.

.

 #ma [op=tom; le=; siz=; n=; b=;]

(tomography). Эта операция выполняет преобразование матрицы в матрицу, используемое в томографии. Исходная матрица имеет размерность [le]*[siz]. Она должна быть записана по столбцам в реальном массиве, начиная с индекса [n]. Порядок записи 11 21 31 ... 12 22 ... Вычисляемая матрица квадратная и имеет размерность [le]*[le] и будет записана в массиве r(), начиная с индекса [b]. Для получения правильного результата необходимо, чтобы матрицы не перекрывались. Элементы новой матрицы с координатами xn,yn отсчитываемыми от середины области изменения индексов, получаются из элементов старой матрицы следующим образом. Каждый столбец старой матрицы ставится в соответствие значению угла A=PI*m/[siz], m=0,....[siz]-1. Координата каждого столбца получается из координат xn,yn по закону x=xn*cos(A)+yn*sin(A). Для этой координаты определяются два ближайших элемента, между которыми производится интерполяция. Если индекс элемента выходит за рамки, то элемент равен нулю. Полученные таким образом значения для всех столбцов суммируются.

.

 #ma [op=gis; b=; le=; n=; sa=;]

(histogram). Эта операция вычисляет гистограмму картинки, записанной как числовой массив. Числовой массив находится в реальном массиве r() начиная с элемента, номер которого [b], а число элементов [le]. Массив должен быть заранее специально приготовлен так, чтобы его элементами были только целые числа от нуля до [n]-1. В принципе, это не обязательно, программа будет работать в любом случае, но ответ будет неправильный. Однако вполне достаточно, чтобы элементы массива были положительными и не больше, чем [n]-1. Программа просматривает весь массив, у каждого элемента берет целое значение снизу и суммирует все такие значения в элементе нового массива, индес которого равен этому целому значению. Исходно все элементы этого массива обнуляются. Результат записывается снова в реальный массив r() начиная с элемента, номер которого [sa]. Размер нового массива [n].

.

 #ma [op=box; wid=; hei=; sty=; bot=; top=; uni=; b=;]

(box). Она выполняет вспомогательные вычисления координат точек для рисования осей прямоугольного графика, внутри которого будет нарисована функция. Такие вычисления проводятся при выполнении операций научной графики в командах #w и #g . Эти операции будут описаны в последующих главах. Данная операция позволяет написать программу, аналогичную #w [op=pf;], но в рамках универсальной графики #eg с возможностью комбинировать много графиков и других элементов непосредственно на ACL, что (хотя и в сложнее для программиста) открывает большие возможности создания графиков типографского качества. Перед вызовом данной операции должны быть заданы параметры [wid], [hei], определяющие область графика (осей), [sty], [bot], [top], задающие видности осей и размеров рисок, аналогично указанным выше командам, [uni] - может быть нулем для автоматической разметки осей и не нулем для ручной разметки, [b] - указывает первый индекс массива r(), в котором должны быть заданы вещественные значения. Если [uni]=0, то таких значений 4, а именно, xmi, xma, fmi,fma, а если [uni]=1, то таких значений 11 и они точно те же, как в команде #w [op=pf;] и других, где размечаются оси. В результате выполнения операции i(1) равно полному числу точек для рисования линий, а i(2) равно полному числу текстов с числовыми значениями при метках. Обозначим эти числа временно как N и M. Тогда числа r(1), ..., r(2*N) содержат координаты точек для линий по правилу x1,y1,x2,y2,.... Числа i(3), ..., i(N+1) указывают на видность сегментов. Если i(3)=1,, то первый сегмент (между 1-й и 2-й точками) соединяется линией, если i(3)=0, то не соединяется и сегмент невидим. Если i(4)=1,, то второй сегмент (между 2-й и 3-й точками) соединяется линией и так далее. Далее в числах r(2*N+1), ..., r(2*N+4*M) содержится информация о текстах. На каждый текст 4 числа. Первые два - это координаты, в которых текст должен быть помещен, третье - это число, которое должно быть показано текстом, четвертое - указание к какой оси текст относится (0 - левая ось, 1 - нижняя, 2 - правая, 3 - верхняя). В зависимости от принадлежности к оси текст ставится либо началом, либо концом, либо серединой в указанную точку. Наконец, 4 числа начиная с r(2*N+4*M+1) указывают xbeg, cx, ybeg, cy - начальные значения осей и масштабные коэффициенты. Эти данные используются для рисования кривой функции. Примером применения этой операции может служить суперкоманда ##pfFCH.

Другой тип специальных вычислений связан с моделированием движения твердых шаров внутри заданного объема. Эту задачу можно полностью решить на ACL, но если шаров много, то скорость вычислений неудовлетворительная. Дело в том, что такой процесс реализуется с использованием команды анимационной графики и в перерыве между быстрой сменой кадров желательно выполнить все вычисления быстро, то есть необходима максимальная скорость. Поэтому написаны три небольшие программы на Java, которые запускаются в работу как отдельные операции. Предполагается, что число шаров задано в переменной N, а диаметр шаров -- в переменной D. Соответственно в первых 4*N элементах реального массива r() записаны параметры для каждого шара, а именно, две координаты x,y центра и две компоненты X,Y скорости (точнее приращение координат за одни шаг). Первая операция

.

 #ma [op=bsm;]

(ball system movement). Она просто выполняет передвижение шаров на один шаг, то есть для каждого шара его координаты меняются по закону x=x+X, y=y+Y. Вторая операция

.

 #ma [op=bsi;]

(ball system interaction). Она выполняет расчет возможного взаимного соударения системы шаров. Алгоритм вычислений стандартный - две частицы, для которых расстояние между их центрами меньше диаметра, обмениваются проекциями скоростей на ось, соединяющую их центры. Проверяются все комбинации пар частиц. Использование этой операции заметно ускоряет процесс. Третья операция

.

 #ma [op=bsr;]

(ball system reflection). Она выполняет расчет возможного отражения шаров от границ области, в которой они находятся. Предполагается, что область описывается ломаной линией, соединяющей точки на плоскости. Эти точки x,y необходимо задать в реальном массиве r() начиная с индекса, который является значением переменной K. Число прямолинейных участков должно быть записано в переменной N1, а число заданных точек должно быть на 1 больше. Если область имеет достаточно сложный контур, то использование этой операции также существенно ускоряет процесс.

      СПЕЦИАЛЬНЫЕ ВЫЧИСЛЕНИЯ В ВИДЕ ПРОЦЕДУР

В языке ACL есть дополнительный набор готовых вычислений, которые не оформлены как операции команды #ma , а вызываются как обычные процедуры языка, но этим процедурам не надо писать куски ACL кода в скобках #pro ... @ . Они имеют специальные названия, которые интерпретатор понимает сразу и сразу выполняет. Код этих процедур написан на языке Java, аналогично всем процедурам, которые вызываются при исполнении команд ACL. Признаться, я даже не могу объяснить почему я не сделал специальную команду, а решил вызывать такие процедуры через команду #e . Но так сделано. Причина видимо в том, что указанные программы являются частью моей личной научной работы и они не являются элементами языка ACL для всех.

По этой причине я совсем не собирался упоминать про них в описании ACL. Многие из указанных программ достаточно сложны и имеют собственное описание. Однако среди них есть такие, которые могут пригодиться кому-то еще кроме меня. По этой причине я решил частично дополнить описание и рассказать про синтаксис таким программ. Эти процедуры имеют название, которое начинается с символа \ (юникод 92), после которого идет номер из двух разрядов (ноль писать обязательно), вот пример

 #e [определение параметров] _\nn

Параметры в виде реальных чисел передаются в программу либо через элементы массива r(), либо через переменные, а целые числа могут передаваться либо через целый массив, либо через параметры. Некоторые программы имеют доступ к переменным и массивам языка, а другие вообще написаны независимо и обмен информацией с ними проводится через запись и считывание файлов. Каждый номер вызывает конкретную программу из тех, что уже написаны. Программы решают разные задачи, которые мне приходилось решать в своей практике. Описание этих задач нужно читать в другом месте. Некоторые из программ настолько сложны, что их описание занимает объем целой статьи. Кроме того, для их использования надо хорошо знать раздел физики, называемый рентгеновская оптика. Тем не менее, я для себя буду здесь записывать короткий справочник по использованию процедур.

\00 -- расчет коэффициента прохождения для многокристального монохроматора на плоскости энергия-угол
\01 -- программа рентгеновской инлайн оптики, т.е. перенос волновой функции через объекты и воздух
\02 -- программа метода стоячих волн, угловая зависимость рентгеновского отражения и выхода вторичных процессов
\03 -- программа расчета лауэ-дифракции в кристалле методами геометрической оптики
\04 -- расчет набора стандартных функций для использования в ACL программах одномерной инлайн оптики
\05 -- программа трехволновой дифракции в кристаллах
\06 -- программа четырехволновой дифракции в кристаллах
\07 -- программа расчета реальных функций двух аргументов (матриц)
\08 -- программа расчета произвольных функций
\09 -- программа spexpro подключения готовых системных программ
\10 -- программа extpro произвольного расчета

Программа 08 является заменой программам 04 и 07, но их пришлось оставить для совместимости со старыми версиями ACL программ. Программы 09 и 10 сделаны специальным образом. Первоначально они задумывались как самостоятельные, без связи с переменными и массивами ACL, то есть они должны были получать входные данные из файлов, а результат снова записывать в файлы. Например, программа 10 может записывать входные данные в файл in.txt в папке tmp, а результат записывать в файл out.dat в той же папке. Это позволяло разрабатывать их в отдельной программе и компилировать отдельно. Затем готовые файлы spexpro.class и extpro.class было достаточно поставить в папки com/vks и com/vku соответственно, и программа уже может работать как внутренняя программа ACL без перекомпиляции всего кода. При этом может существовать много версий файлов spexpro.class и extpro.class, созданных в отдельных программах.

Через некоторое время я сообразил, что программы можно также разрабатывать отдельно симулируя те же массивы, которые используются в ACL. При этом после подключения к ACL они могут использовать память языка, что более удобно чем файлы. Так что сейчас используются оба способа. Кроме того, сейчас сделано так, что эти программы имеют аргумент в виде целого числа, совпадающего с параметром [n] при вызове, то есть их надо вызывать так
#e [n=1;] \09
А сами специальные программы spexpro и extpro просто вызывают конкретную программу из списка программ согласно значению параметра [n]. Отдельные программы надо писать таким образом, чтобы код spexpro.java и extpro.java находился в папках com/vks и com/vku соответственно, а в остальном они могут быть произвольными. Фактически такой способ является универсальным и позволяет заменить все другие. Надо было сразу так сделать. Но неразумно переписывать старый код, поэтому этот способ ориентирован на будущее. Например, по номеру [n=1;] \09 запускается мой большой текстовый редактор, который раньше был написан как отдельная программа. И я, фактически только так его и использую последнее время. А у программы \10 уже много значений параметра [n] используется.