text:
# Angular Resolution of the Pachmarhi Array of Čerenkov Telescopes

## 1 Introduction

The $`\gamma `$ray sky has been successfully observed by satellite based missions in the energy range between 50 $`keV`$ and 30 $`GeV`$ . At higher energies the steeply falling $`\gamma `$ray spectrum makes it impossible to observe by satellite based detectors, primarily because of the limited detector size and exposure time available for a given source. Alternatively, the ground based atmospheric Čerenkov technique has been successfully exploited at these energies. Over the years, this technique has proved to be the most sensitive in studying the celestial $`\gamma `$rays in the very high energy (VHE) range, $`100GeV`$ \- 50 $`TeV`$.

Celestial VHE $`\gamma `$rays initiate an electromagnetic cascade in the atmosphere as they enter the earth’s atmosphere. The electrons and positrons in the cascade, being relativistic, emit Čerenkov light as they propagate down the atmosphere resulting in a faint flash of bluish light lasting a few nanoseconds. This directional Čerenkov flash is detected on the ground by conventional optical light detectors during moon-less clear nights. Taking advantage of the fact that the Čerenkov flash is highly collimated, within a cone of half angle $`1^{}`$, along the direction of the incident particle most experiments simply limit the optical field of view of the Čerenkov telescopes to a small region of the sky.

However the main drawback of this technique is the presence of abundant cosmic ray background which severely limit the sensitivity of this technique. Numerous efforts to develop methods to distinguish the Čerenkov light pool produced by cosmic $`\gamma `$rays from that by the cosmic rays led to two important techniques based on complementary ways of viewing the cascade, viz. the angular sampling or the imaging technique and spatial sampling or the wavefront sampling technique. Both these techniques are currently being employed by different groups .

The possibility of timing the Čerenkov light front for the estimation of arrival direction of the primary has been realized as early as 1968 by Tornabene and Cusimano . It was later demonstrated by the Durham group  that while viewing Crab pulsar the on-axis events were richer in $`\gamma `$-rays than off-axis events. It was also used by Gupta et al. to improve the signal to noise ratio of atmospheric Čerenkov telescopes consisting of large mirrors with poor optical quality. The basic principle of all these experiments is the fast timing technique using spatially separated Čerenkov telescopes.

The signal to noise ratio in such an experiment is given by

$$\frac{S}{N}\sqrt{\frac{AT}{(1f_r)F_c\omega }}f_\gamma F_\gamma $$
(1)
where $`F_c`$ and $`F_\gamma `$ are the fluxes of omni-directional cosmic rays and $`\gamma `$ rays from a point source respectively. A is the collection area of the array, T is the time of observations, $`\omega `$ is the solid angle of acceptance, $`f_r`$ is the fraction of showers due to cosmic rays that are identified and rejected as background and $`f_\gamma `$ is the fraction of showers due to $`\gamma `$rays that are identified as signal and hence retained. In order to achieve high $`\frac{S}{N}`$ apart from increasing $`f_\gamma `$ and $`f_r`$ one could either increase the collection area and the observation time or decrease the solid angle of acceptance. For a given exposure time and available hardware resources one can possibly increase S/N by only reducing $`\omega `$ as celestial $`\gamma `$ rays from point sources are directional while cosmic rays are isotropic. Due to the finite opening angle of the Čerenkov cone and the spread in the arrival angle of Čerenkov photons the aperture of the telescopes have to be restricted to few degrees, which sets a lower limit to $`\omega `$. However, it is possible to improve the $`\frac{S}{N}`$ for point sources without losing Čerenkov light if the direction of arrival of primary particles is estimated accurately . The shower axis retains the original direction of the primary.

To estimate the arrival direction accurately, the error in the reconstructed direction, i.e. the angular resolution has to be very small. The two dominant factors which contribute to the angular resolution are the average distance, D, between the telescopes and $`\delta `$t, the uncertainty in the measurement of arrival time of photons at the telescopes. For example, the angular resolution in the zenith angle, $`\theta `$ for $`n`$ detectors is given by :

$$\delta \theta =\frac{c\delta t}{D\mathrm{cos}\theta }\sqrt{\frac{2}{n}}$$
(2)
Therefore, a large number of telescopes with instrumentation to measure the relative arrival time of photons and separated by large distances are needed to reconstruct the shower front and estimate the direction of arrival of the shower fairly accurately. On the other hand, in case of the imaging technique the angular resolution is limited by PMT(pixel) size which is typically of the order of a quarter degree. The improvement in the signal to noise ratio by restricting the angle of acceptance to $`\delta \theta `$ would be:

$$\sqrt{\frac{\omega }{\pi \delta \theta ^2}}$$

This is a very significant advantage of non-imaging arrays with high angular resolution. Moreover, this method of reducing background is independent of the primary species. Hence it is extremely useful at very low primary energies where the cosmic ray electrons form a significant source of background. Electrons cannot be discriminated easily, unlike hadrons, since they too undergo electromagnetic interactions in the atmosphere akin to $`\gamma `$ rays. However one has to bear in mind that the angular resolution could be poorer at lower primary energies due to paucity of Čerenkov photons at the ground level.

In the rest of the paper we discuss the method of obtaining the arrival direction of the primary radiation using our array of Čerenkov telescopes and estimate its error.

## 2 Pachmarhi Array of Čerenkov Telescopes

The experiment at Pachmarhi (longitude: 78 26 E, latitude: 22 28$`{}_{}{}^{}N`$ and altitude: 1075 $`m`$), is based on the wavefront sampling technique and employs an array of Čerenkov telescopes, called the Pachmarhi Array of Čerenkov Telescopes(PACT), to sample the Čerenkov light pool. PACT is now fully functional.

The experimental set-up of PACT has been explained in detail elsewhere . Briefly, it consists of a $`5\times 5`$ array of atmospheric Čerenkov telescopes deployed in the form of a rectangular matrix with a separation of 25 m in the N-S direction and 20 m in the E-W direction. Each telescope consists of 7 parabolic mirrors of 0.9 m diameter mounted paraxially and having a focal length of 90 cm. Each mirror is viewed by a fast photo-multiplier tube (PMT, EMI 9807B) at the focus behind a circular mask of $``$ 3 diameter. However the field of view is limited by the diameter of the PMT photo-cathode to $``$ 2.9 FWHM.

The array has been divided into 4 sectors with six telescopes in each. The pulses from 7 PMTs in a telescope are added linearly to form a telescope sum pulse called royal sum. Each royal sum from the 6 telescopes in a sector are suitably discriminated (typical royal sum rates $``$ 30-50 kHz.) and a trigger is generated by a coincidence of any 4 of these 6 royal sums. The typical event rate is $``$ 2-5 Hz. For every event, information on the relative arrival times and density of Čerenkov photons are recorded for the 6 peripheral mirrors/PMT in each telescope in each sector. The relative arrival times of royal sum pulses are recorded both in the respective sector and in the central data processing station.

Thus, PACT measures the arrival time of shower front at various locations within the Čerenkov light pool at two distance scales, short range (intra-telescope) and long range (inter-telescope). The arrival direction of the shower is estimated from the measured arrival times of Čerenkov photon front at each of the spatially separated telescopes while the six adjacent measurements in a given telescope could be used to study the fluctuations in the measured shower parameters. The dispersion of photon arrival times at each telescope contain the signature of the primary species  and hence could be used to distinguish between $`\gamma `$ rays from the background . The density measurements on the other hand enable us to estimate the energy of the primary species as well as to reject cosmic ray background .

## 3 Alignment of Mirrors and Telescope

### 3.1 Pointing Accuracy of Telescopes

The telescopes are equatorially mounted and each telescope is independently steerable in right ascension and declination within $`\pm 45^{}`$. The movement of the telescopes is remotely controlled by Automatic Computerized Telescope Orientation System (ACTOS) . The hardware for ACTOS consists of a semi-intelligent closed loop stepper motor drive system which senses the angular position using a gravity based angle transducer called clinometer with an accuracy of 1. The two clinometers, one each in N-S and E-W direction, are accurately calibrated using stars at various hour-angles and declinations.

In order to estimate the source pointing error of the system, the telescopes were oriented to different bright stars at random positions in the sky and the offsets in orthogonal directions (ascension and declination) with respect to the star were noted. Using this data it is concluded that the system can orient the telescopes to the putative source with a mean offset of $`0^{}.003\pm 0^{}.1`$. The source pointing is monitored with an accuracy of $`0^{}.05`$ and corrected in real time during tracking. The uncertainty in pointing represents the subjective error in manual reading of the star position at the cross wire of the guiding telescope field of view ($`\pm 0^{}.05`$).

### 3.2 Measurement of Alignment Accuracy of Mirrors

As the seven mirrors of a telescope are mounted on a single mount it is necessary to ensure that all their optic axes are parallel to each other so that they view the same part of the sky. This alignment of mirrors is done manually by positioning the star image at the centre of the focal plane.

To check the accuracy of alignment of the mirrors and telescopes a bright star drift scan is carried out. The telescope is aligned to an isolated bright star (typically of visual magnitude 2 to 3). Then the telescope is moved to the west by $`4^{}`$ and at a suitable time the telescope tracking is switched off. The counting rates from each of the PMTs are monitored every second and recorded. The count rates stay reasonably constant until the star walks into the field of view when they increase, go through a maximum and return to the background value as shown in figure 1. The background count rates before and after the star transit are fitted to a linear function. The background is subtracted from the count rates during the star transit by interpolation. The background corrected count rates are then fitted to a quadratic function. Figure 1 shows one such fit to a typical count rate profile of a mirror. Table 1 shows the summary of results of a typical bright star scan. It shows in column 3 the FWHM of the drift scan profile of the count rates due to the star. The offset of the centroid of this profile with respect to the centre of the field of view (i.e. expected transit time) is shown in the last column. This offset is the error in the alignment in right ascension. The relative FWHM’s of the profiles could indicate the misalignment in declination, if any. On the other hand the absolute values of FWHM of the count rate profiles are a function of the PMT gains, image quality (point source image size $`1^{}`$), star brightness $`etc`$. Using this method, it is ensured that the optical axes of all the 7 mirrors in a telescope are parallel to each other within an error of $`0^{}.2`$ or less. If the error exceeds this value the particular mirror is re-aligned and re-checked. This method of alignment is similar to that adapted for aligning the imaging telescope arrays .

## 4 Estimation of Timing Resolution

In order to estimate the angular resolution one has to estimate the error in timing measurement. To determine the arrival time of photons accurately, the PMT should have a fast response time and low dispersion in the photo-electron transit times. We use a low noise, high gain tube having rise-time of $``$ 2 ns, pulse width of 3 ns and transit time jitter of 1.9 ns (FWHM). The PMT signals are brought to the recording stations through low loss coaxial cables (RG 213). The intrinsic timing jitter of the signals from the PMTs limit the accuracy of timing measurements.

The timing resolution is estimated from the data. The relative delays of 6 peripheral PMT’s of each telescope and royal sum pulses with respect to the trigger pulse were recorded using the fast time to digital converters(TDC). The signals were delayed suitably using ECL based delay generators having an intrinsic stability of $`\pm 0.1ns`$. Linearity of TDC was monitored by periodically calibrating them with standard cable delays. The LeCroy(Model # 2228A) and Philips Scientific(Model # 7186) TDC modules that are used in the experiment are set to a sensitivity of 0.25 and 0.2 ns per count, respectively.

Data are collected with all telescopes pointing to zenith. The width of the distribution of differences in arrival times of signals ($`\delta t_{ij}`$) of two TDC channels is an indication of the limiting accuracy of timing measurement, provided the signals originate from PMT’s located nearby. The distribution of $`\delta t_{ij}`$ tends to broaden due to one or more reasons mentioned below:

(a) fluctuations in the arrival time and density of Čerenkov photons.

(b) errors in mirror pointing result in sampling different angular regions, which in turn leads to pulse slewing.

(c) spatial separation between the mirrors/telescopes. The mean arrival time as well as its fluctuations is a function of core distance .

(d) fluctuations in electron transit time in PMT’s, the differences in the cable delays and the propagation delays in the electronic circuits.

Delays due to unequal cable lengths, and propagation delays in the electronic circuits have been matched to an accuracy of about 0.25 ns. Mirrors are aligned to within $`\pm 0^{}.2`$ as mentioned earlier. The effect of slewing responses of the PMT’s is often determined by using filter wheels to attenuate laser calibration light flashes transported to each mirror through optical fibers . No such measurements were made in this case since the pulse rise time is still reasonably small ($`2.4ns`$). The effects due to finite separation between the detectors is corrected as discussed below and shown in figure 2.

In the case of signals from individual PMT’s of a telescope the distances between the detectors is of the order of a metre. Only those combinations corresponding to neighboring PMT’s are considered in order to minimize the contributions from fluctuations in the arrival time of Čerenkov photons and distance dependent effects, to the uncertainty in timing measurement. The results are shown in the table 2. The limiting accuracy of timing measurement is about 1 $`ns`$ for individual PMT signals.

In the case of royal sum pulses, the distance between neighboring telescopes is 20 m or more, which is too large to be ignored. Moreover, the spread in the arrival time depends upon how well the signals of a telescope are added. To remove the distance dependent effects, the standard deviation $`\sigma _{ij}`$ are plotted as a function of separation between the telescopes, $`i`$ and $`j`$ and is shown in figure 2. In the past, the latter effect was minimized by tilting the telescopes towards each other by about a degree and the timing resolution was measured . However we feel that the present method of reducing contribution from the detector separation is more accurate. From the figure it can be seen that the error $`\sigma _{ij}`$ clearly increases with the separation as expected. The FWHM for two telescopes side by side is estimated by extrapolation, thus removing the effect of finite separation between the telescopes on the $`\sigma _{ij}`$ distribution. The corresponding values are shown in table 2. These values are comparable to the intrinsic timing jitter for proton primaries estimated from simulation studies for PACT telescopes . This means that the contribution to timing resolution from instrumental effects is negligibly small.

## 5 Estimation of Arrival Direction of Shower

The arrival direction of a shower is determined by reconstructing the shower front using the relative arrival time of Čerenkov photons at each telescopes. We fit the Čerenkov photon front to a plane using the measured arrival times, the normal to this plane gives the direction of shower axis.

### 5.1 Calculation of Time-offsets

The relative arrival time of pulses as measured in the experiment is not the relative arrival time of Čerenkov photons at the PMT, needed for reconstructing the shower front. A finite but constant delay between pulses from different PMT’s arise due to unequal cable lengths, differences in electronic propagation delays and differences in photo-multiplier transit time etc. These are termed as T0 or time-offsets. Thus the measured relative arrival times have to be corrected for this time-offsets to get the relative arrival time of Čerenkov photons at the PMT. In the absence of a calibration light source which could generate pulses in all PMT’s simultaneously, these time-offsets have been determined from the data itself .

The time difference between Čerenkov photons arriving at PMT’s (or telescopes) located side by side should be zero on the average for vertical showers. Thus the average of time delays between two PMT’s, from a large sample of data, is entirely due to difference between the two time-offsets. We chose vertical showers in order to eliminate the errors in the estimation of offsets due to finite telescope separations for inclined showers.

If $`T0_i`$ and $`T0_j`$ are the time offsets for the PMTs i and j, we can write an equation of the form

$$(T0_iT0_j)=C_{ij}$$
(3)
where $`C_{ij}`$ is the mean delay between a pair of PMT’s i and j after correcting for the time difference due to differences in height (z-coordinates) of PMT’s(or telescopes) .

The total $`\chi ^2`$ can be expressed as

$$\chi ^2=\underset{i,j=1;ij}{\overset{n}{}}W_{ij}(T0_iT0_jC_{ij})^2$$
(4)
where $`n`$ is the total number of PMT’s, $`W_{ij}`$’s are the statistical weight factors, $`(W_{ij}=1/\sigma _{}^{C}{}_{ij}{}^{}{}_{}{}^{2}`$, where $`\sigma _{}^{C}{}_{ij}{}^{}`$ is the uncertainty in determining $`C_{ij}`$). $`\chi ^2`$ minimization will give a set of $`n`$ equations of the form

$$\underset{j=1;ij}{\overset{n}{}}W_{ij}(T0_iT0_j)=\underset{j=1;ij}{\overset{n}{}}W_{ij}C_{ij}=C_i^{}$$
(5)
Thus we have a set of simultaneous equations involving W, T0 and $`C^{}`$. All these sets of equations can be written in the form of matrix equation

$$WT0=C^{}$$
(6)
Solving this matrix equation one gets an estimate of time-offsets for (n-1) detectors with respect to $`n^{th}`$ detector. It is to be noted that the differences in time delays due to the differences in the heights of the telescopes were taken off while calculating these time-offsets, as already mentioned.

However the mean time difference $`<\delta t>`$ (ns) is a function of core distance R (m) and varies, for vertical showers, as

$$<\delta t>=\left(\frac{\delta d}{Hc}\right)R$$

where $`\delta d`$ is the separation between the telescopes, H the radius of curvature of the Čerenkov shower front and c the velocity of light. Typically, $`<\delta t>=0.0117R\mathrm{𝑛𝑠}`$ for the energy and altitude concerning us (H $``$ 10 km). In our experiment, the average separation between the telescopes in a sector is 35 m-50 m whereas the separation between the PMT’s in a telescope is of the order of a metre. This can introduce a systematic error in the estimated values of T. For a proton shower of average energy ($`1.8TeV`$) the mean collection radius is $`100m`$. Using this it is estimated that the maximum systematic error in T arising from the Čerenkov shower front curvature is $`1.2ns`$.

### 5.2 Reconstruction of Arrival Direction

Using the plane front approximation the arrival direction of the shower is estimated as follows . If $`x_i`$,$`y_i`$,$`z_i`$ are the coordinates of the $`i^{th}`$ PMT, $`l,m,n`$ the direction cosines of the shower axis and $`t_i`$ the arrival time of the photons at this PMT then the equation relating them is given by

$$lx_i+my_i+nz_i+c(t_it_0)=0$$
(7)
where $`t_0`$ is the time at which the shower front passes through the origin of the coordinate system. Then the arrival direction of the shower can be estimated by a $`\chi ^2`$ minimization where

$$\chi ^2=\underset{i=1}{\overset{n}{}}w_i(lx_i+my_i+nz_i+c(t_it_0))^2$$
(8)
where $`w_i`$’s are the statistical weight factors for the $`i^{th}`$ timing measurement(= 1/$`\sigma _i^2`$ where $`\sigma _i`$ is the uncertainty in the timing measurement for the $`i^{th}`$ detector). The values of $`l,m,n`$ and $`t_0`$ are calculated using the equations $`\chi ^2/l=0`$, $`\chi ^2/m=0`$, $`\chi ^2/t_0=0`$ and $`l^2+m^2+n^2=1`$. We have ignored terms containing $`n/l`$ and $`n/m`$ as these are small. The contribution to the error in angle determination is negligibly small up to an angle of $`80^{}`$ under this assumption.

In the first iteration all PMT’s with valid delays are used for the fit and the observed - expected delays are calculated for each. In the subsequent iteration, the PMT with the largest deviation is rejected if the absolute difference between observed and expected delays is greater<sup>1</sup><sup>1</sup>1$`\pm 3\sigma `$, The standard deviation of observed - expected delay $`1ns`$  than 3 ns. This is to minimize the effect of isolated large fluctuations in the photon arrival times . Figure 3 shows the observed - expected delay distributions for various telescopes, based on royal sum pulses. This process of fitting is continued until there are no more large deviations or the number of available PMT’s have become less than 4, in which case the shower front reconstruction process is aborted.

Figure 4 shows the zenith and azimuthal angle distributions of the reconstructed arrival directions using the procedure explained above. The reconstruction of the shower front was carried out using 25 telescope data collected with telescopes in the vertical position. The observed distributions are consistent with those expected.

Further, it has to be ascertained that the telescopes are pointing to the true zenith direction. A tilt in the observer’s co-ordinate system will imply a systematic error in pointing to the source direction. It is ensured that the assumed horizontal plane is really perpendicular to the direction of true zenith from the projected angle distributions of the fitted directions. The projected angles $`\alpha `$ and $`\beta `$ are related to the zenith angle $`\theta `$ and azimuth $`\varphi `$ by

$$tan\alpha =tan\theta cos\varphi $$

and

$$tan\beta =tan\theta sin\varphi $$

where $`\alpha `$ and $`\beta `$ are the projected angles in the N-S and E-W plane respectively. The frequency distributions of $`\alpha `$ and $`\beta `$ are shown in figure 5. The ratio of the number of events above and below zero are 1.0 $`\pm `$ 0.03 and 0.96 $`\pm `$ 0.03 for $`\alpha `$ and $`\beta `$ respectively. The symmetry of these distributions thus verify the orthogonality of the telescope pointing direction to the horizontal.

## 6 Angular Resolution of PACT

We define the angular resolution of the array by the angle which corresponds to the 68% acceptance of mainly proton events. This definition is consistent with that for the HESS array except that their definition refers to $`\gamma `$-ray showers only . The angular resolution of PACT has been estimated by using the split array method. The array is divided into two independent parts and the arrival direction is estimated for each shower from these two sub-arrays. The distribution of space angle ($`\psi `$) between these two estimates is a measure of the angular resolution, i.e. the accuracy with which one can estimate the arrival direction. Since there are two independent estimates of the direction, the angular resolution is given by the peak<sup>2</sup><sup>2</sup>2for small angle approximation $`\sigma _\psi `$ of the point spread function coincides with the peak of the space angle ($`\psi `$) distribution of the distribution of space angle between the two directions as $`\frac{Peak}{\sqrt{2}}`$.

Similarly, the uncertainties in the zenith angle $`\theta `$ and azimuth $`\varphi `$ are estimated from the distributions of $`\mathrm{\Delta }\theta `$ and $`\mathrm{\Delta }\varphi `$, the differences between the corresponding angles from two independent estimates for a shower. Figure 6 shows these distributions for $`\mathrm{\Delta }\theta `$ and $`\mathrm{sin}\theta .\mathrm{\Delta }\varphi `$, respectively.

### 6.1  Using Royal Sum pulses

The array is divided into two independent parts of equal number of telescopes, say 12 telescopes from sectors 1,2 and sectors 3,4. The arrival direction is estimated for each shower from these two independent arrays using the relative arrival times of royal sum pulses. The distribution of the space angle between the two estimated directions is shown in figure 7. Several combinations were tried to understand the dependence of angular resolution on the separation between the telescopes and/or the number of detectors. The results are summarized in table 3. The term ‘odd-even’ refers to 3 telescopes each from Sector 3 and 4 grouped into one set of 6 telescopes and the remaining 6 grouped into another.

It is seen from table 3 that the angular resolution improves as the separation between the detectors increases and also as the number of detectors increases, as expected. The dependence of angular resolution on D, the distance between telescopes and n, the number of telescopes used in the fit are shown in figure 8. The panel (a) of the figure shows the variation of angular resolution with D for (i) 4 telescopes and (ii) 12 telescopes ($``$ and $``$ respectively). The dependence on D is calculated to be $`\sigma _\psi D^{0.94\pm 0.08}`$ and $`\sigma _\psi D^{1.05\pm 0.05}`$ respectively, which is compliant with theoretically expected 1/D dependence.

Once the 1/D dependence is established, the angular resolution is then corrected for the distance to extract the dependence on the number of degrees of freedom. This is shown in the panel (b) of figure 8. It is seen that $`\sigma _\psi n^{0.41\pm 0.11}`$, consistent with what one expects, on a simple statistical basis, that $`\sigma _\psi `$ will decrease as $`\frac{1}{\sqrt{n}}`$.

It is seen from table 3 that the best value of angular resolution, obtained using 12 royal sum pulse information, is $`0^{}.33`$. Since we have 25 telescopes, there will be a further improvement in the angular resolution by a factor of $`\sqrt{2}`$ when the arrival direction is estimated using all the telescopes in the array. Thus the angular resolution of the array is estimated to be $`0^{}.24`$, from the royal sum information of all telescopes. The corresponding threshold energy of the primary $`\gamma `$ray is about 1.5 $`TeV`$.

### 6.2 Using Individual PMT timing pulses

The information on the relative arrival times of Čerenkov photons at the individual PMT’s are available only within a sector. Only those events with information in all 6 telescopes in sectors 3 and/or 4 are used for estimating angular resolution. It is obtained by dividing the data from the sectors into 2 subsets and obtaining the space angle between the two arrival direction estimates from the two subsets of data. The distributions of deviations of zenith and azimuth angle estimates are shown in figure 9 which corresponds to data from sector 3.

Information for the 3 PMT’s(labelled A, C and E) are grouped into one set while those for the remaining 3 PMT’s of a telescope(B, D and F) are grouped into another set. The results are summarized in table 4 . The rows 1 and 4 correspond to the cases in which two subsets are obtained by demanding all A PMT’s of 6 telescopes, all B PMT’s etc. Rows 2 and 5 correspond to cases in which the 1st group consists of at least one valid TDC information in any of 3 A, C, E PMT’ and the second group from any one of 3 B, D, F PMT’s. Similarly the rows 3 and 6 correspond to two sets with at least 2 valid TDC’s in each telescope. Finally, row 7 refers to the case in which the arrival directions are obtained separately from sector 3 and sector 4 events and collating event arrival times to pick common events. Column 2 shows the corresponding number of detectors used in the fit for all cases. Figure 10 shows the distribution of space angle between the two direction estimates made using individual PMT signals in a sector.

A better estimate of arrival direction of the shower is possible with the use of individual PMT signals as seen in table 4. The angular resolution obtained using individual PMT’s of just one sector is $`0^{}.24`$ (row 6 of table 4). When one considers the entire array the average separation doubles while the number of degrees of freedom quadruples. Consequently, there will be an improvement by a factor 4 in the accuracy of the arrival angle estimate. Based on these considerations, the angular resolution of the entire array is estimated to be, $`\sigma _\psi `$ $``$ $`0^{}.06.`$, for $``$ 2 active PMT’s per telescope<sup>3</sup><sup>3</sup>3It may be noted that the average number of timing signals from individual mirrors in a sector is 15 when we demand $``$ 2 PMT’s and 27 when we demand $``$ 4 PMT’s in each telescope..

We have 6 timing information from each telescope. If we use the timing information from a totality of 150 PMT’s in the entire array, then the angular resolution works out to be:

$$\sigma _\psi 0^{}.04$$

The same figure is arrived at when we use row 7 of table 4. This estimate is based on the fact that $``$ 4 PMTs have valid TDC information in a telescope. The primary energy threshold for such events is relatively higher, about 3 $`TeV`$. The errors in the arrival direction estimate are larger for lower energy showers. While the angular resolution of the imaging telescopes are limited by the PMT sizes, only the future imaging telescope arrays (like VERITAS  or HESS ) could achieve a better angular resolution. PACT is able to achieve this because it is a distributed array with multiple sampling of Čerenkov light pool.

## 7 Discussions

For a wavefront sampling experiments such as PACT, the technique by which off-axis events are rejected efficiently using their fine angular resolution, is vital. Hence significant effort has gone in to the improvement of the angular resolution of such arrays. In the case of the imaging telescopes the angular resolution is mainly decided by the size of the photo-tubes in the imaging camera at the focal point. However a new technique of stereoscopic viewing has been developed in 1996 by HEGRA group wherein the same shower is imaged by multiple telescopes. Using this technique the shower axis can be completely reconstructed in space using the information concerning the location of the images within each cameras and their angular orientation . The reconstruction procedure allows determination of the shower direction, the core location and the height of the shower maximum. This technique has greatly improved the angular resolution of imaging telescope arrays at the cost of reduced collection area. On the other hand the improvement of angular resolution of a wavefront sampling array results in a direct improvement in the flux sensitivity since the already increased collection area of such arrays is unaffected by the angle reconstruction. It may be noted here that a 3-dimensional reconstruction of the shower is possible for the a non-imaging array by making use of the curvature of the Čerenkov light front .

From the results of the present study it is seen that two independent measurements of the angular resolutions are $`0^{}.24`$ and $`0^{}.06`$ resulting from the timing measurements of telescopes and mirrors of the entire array respectively for similar energy events. The essential difference between the two measurements is an increase in the number of degrees of freedom by a factor of $`2.5`$ in the latter case. From this consideration the expected angular resolution for individual mirrors is $`0.24/\sqrt{2.5}=0^{}.15`$ which is almost twice the value of $`0^{}.06`$ mentioned above. These are consistent since the timing uncertainty for royal sums is higher by the same factor.

Figure 11 shows differential plots of the number of events per solid angle $`\frac{dN}{d\omega }`$ as a function of the zenith angle. The panel $`(a)`$ shows plots for events whose arrival directions have been estimated using timing signals from royal sums while the panel $`(b)`$ shows the same for individual mirror timing information. The two curves shown in panel $`(a)`$ correspond to events detected simultaneously in all the 4 sectors (solid line) and in one sector(# 3) only (dashed line). This provides a method of measuring the effective opening angle of the array . From these plots it can be seen that the array field of view is around $`2^{}.3`$. This is significantly lower than the geometric opening angle of $`2^{}.9`$. It may be noted that the effective field of view of the array is expected to be lower than the geometric value primarily because the event trigger efficiency falls with increasing core distance and zenith angle. The other possible factors which could potentially reduce the effective field of view like the residual alignment errors and non-ideal optical quality of the mirrors are relatively less significant. The effective FWHM of the array ($`2^{}.3`$) for which the mean telescope separation is $`70m`$ is lower than that for a sector ($`2^{}.6`$) for which the mean telescope separation is $`32m`$ supporting the above argument.

The data used for estimating the angular resolution of PACT reported here consists of cosmic rays. It is well known that the protons exhibit a larger (by almost a factor of 4) intrinsic timing jitter relative to $`\gamma `$-ray primaries . While the contribution of systematic effects to the timing resolution (which is estimated to be small) is independent of the primary species, the decreased intrinsic timing jitter is expected to improve the angular resolution for $`\gamma `$-ray primaries significantly.

Table 5 summarizes the angular resolution of imaging as well as non-imaging atmospheric Čerenkov telescopes currently in operation and under construction. It can be seen that the angular resolutions of most of the single imaging telescopes are modest, in the range $`0^{}.10^{}.7`$. Wavefront sampling experiments using the modified solar concentrator arrays like the STACEE, CELESTE, GRAAL $`etc.`$, will not be able to exploit their angle reconstruction technique to reject a significant cosmic ray background primarily because of their rather small field of view, in the range $`0^{}.240^{}.7`$ (FWHM) which are comparable to their angular resolutions. PACT, on the other hand, has the best angular resolution among the non-imaging experiments which can be used to enhance the signal to noise ratio significantly. PACT has been able to achieve this because of multiple fast timing measurements at each telescope. This has two-fold advantages: firstly, it provides an increased number of degrees of freedom which in turn improves the accuracy of angle reconstruction and secondly, it provides a means of computing the Čerenkov photon arrival time dispersion at each telescope. The timing jitter is a species sensitive parameter and hence will enable us reject a significant fraction of on-axis hadronic showers. The signal to noise ratio could also be improved by using the normalized $`\chi ^2`$ values resulting from a spherical fit to the timing information of a shower. This is expected to be larger for hadron initiated showers .

The future imaging telescope arrays are expected to achieve unprecedented angular resolution by employing the stereoscopic technique. They are able to reconstruct the shower in 3-dimensions using the multiple images in 3 or more imaging telescopes. From table 5 it can be seen that both the HEGRA and the SHALON-ALATOO arrays are able to improve the angular resolution of their single imaging telescopes by a factor of $`4`$ by using 2 or more of them in stereoscopic mode. It may be mentioned here that the quoted angular resolutions for the future projects listed in the table are from simulation results for $`\gamma `$-ray primaries while the rest are derived from measurements from cosmic ray events.

It is well known that the Čerenkov light front has a curvature. It was seen that when two well separated Čerenkov telescopes were tilted towards each other by about a degree the coincidence rate increased and also reduced the spread in the time separation between them . This indicated, as claimed by the authors, the presence of curvature in the photon front. More recently, it was shown that the radius of curvature of the front is equal to the height of the shower maximum from the observation altitude . Hence it is clear that a plane front approximation of the Čerenkov light front will introduce a systematic error in the arrival angle reconstruction. The large separation of telescopes, to some extent, offsets the worsening of resolution due to the curvature of the light front. It may also be pointed out that the effect of curvature on the angular resolution is more significant for near vertical showers. In the case of inclined showers, the arrival time differences at spatially separated telescopes due to shower axis inclination far exceed the differences arising out of the shower front curvature. However it has been argued that the angular resolution of given array will improve if one corrects for the wavefront curvature . The details of the systematic effects due to the plane front approximation and the improvement in the accuracy of estimated arrival directions when the curvature of the shower front is taken into account are currently under investigation. A paper based on the results of these simulation studies is under preparation.

## 8 Summary

A detailed analysis of the angular resolution of PACT using data collected with telescopes pointing to zenith is presented. The improvement in the angular resolution with larger separation between detectors (D) and with increase in the number of degrees of freedom (n) has been verified. The angular resolution $`\sigma _\psi `$, is found to improve according to the relation:

$$\sigma _\psi \frac{1}{D\sqrt{n}}$$

There are two types of angular resolutions that could be defined for PACT. Firstly, the fast timing information from the 25 telescopes spread in an area of about $`80m\times 100m`$ yield an angular resolution of $`0^{}.2`$. Secondly the timing information from the mirrors constituting the 25 telescopes yield the best angular resolution of $`2.4^{}`$ at a $`\gamma `$-ray energy of around 3 $`TeV`$. This is the best angular resolution achieved for any ground based atmospheric telescope system which will probably be superseded by the proposed VERITAS or HESS imaging telescope array.

The angular resolution for $`\gamma `$-ray primaries could be significantly better than what is presented here, which is based on cosmic ray primaries, because of two main reasons. Firstly, photon arrival time jitter for $`\gamma `$-ray primaries is far less than that for charged cosmic rays and secondly the radius of curvature for $`\gamma `$-ray primaries is more than that for cosmic rays which consequently reduces the effects of the shower front curvature for the former.

We would like to acknowledge the hardwork and dedication of our colleagues who helped us in the installation of PACT, its instrumentation and the subsequent observations: Messrs S. S. Upadhya, K. Gothe, B. L. Venkateshmurthy, B. K. Nagesh, K. K. Rao, A. J. Stanislaus, S. K. Rao, P. V. Sudershanan, P. N. Purohit, S. K. Sharma, M. S. Pose, A. I. D’Souza, S. R. Joshi. We would like to thank the anonymous referee for a careful, patient and thorough job of refereeing this paper.


format: text

text:
# 1. INTRODUCTION

## 1. INTRODUCTION

Real objects in general are fractional-order systems, although in some types of systems the order is very close to an integer order. Such systems are mainly electronic systems composed of quality electronic elements. So far, however, the systems were described as integer-order systems, regardless of the negative consequences caused by neglecting the real order of the system . Disregarding the fractional order of the system was caused mainly by the nonexistence of simple mathematical tools for the description of such systems. Since major advances have been made in this area in the last few years  it is possible to consider also the real order of the dynamical systems. Such models are more adequate for the description of dynamical systems with distributed parameters than integer-order models with concentrated parameters. With regard to this, in the task of identification, it is necessary to consider also the fractional order of the dynamical system. In this contribution we will concentrate mainly on the identification of parameters (including the order of derivatives) for a chosen structure of the model with emphasis on the methods of evaluation of transfer characteristics. The verification of the correctness of parameter identification will be done by using it in systems with known parameters.

## 2. DEFINITION OF THE FRACTIONAL-ORDER SYSTEM MODEL

For the purposes of later development let us consider as the model of a dynamical system the following fractional-order three-member differential equation

$$a_2y^{(\alpha )}(t)+a_1y^{(\beta )}(t)+a_0y(t)=u(t)$$
(1)
where $`\alpha `$ and $`\beta `$ are in general real, with initial conditions $`y^{(\beta )}(0)=0`$ and $`y(0)=0`$.

## 3. NUMERICAL SOLUTION OF THE FRACTIONAL-ORDER
DIFFERENTIAL EQUATION

For the numerical calculation of the fractional-order diffrential equation (1) we use for an approximation of fractional derivatives in the equation (1) the following formula taken from

$$y^{(\alpha )}(t)_{(tL)}D_t^\alpha y(t)=h^\alpha \underset{j=0}{\overset{N(t)}{}}b_jy(tjh),$$
(2)
where $`L`$ is the ”memory length”, $`h`$ is the time step,

$$N(t)=\mathrm{min}\{\left[\frac{t}{h}\right],\left[\frac{L}{h}\right]\},$$

$`[z]`$ is the integer part of $`z`$,

$$b_j=(1)^j\left(\genfrac{}{}{0pt}{}{\alpha }{j}\right)$$
(3)
where $`\left(\genfrac{}{}{0pt}{}{\alpha }{j}\right)`$ is a binomial coefficient.

The approximation of the equation (1) in discrete time steps $`t_m(m=2,3,\mathrm{})`$ has the following form

$$a_2h^\alpha \underset{j=0}{\overset{m}{}}b_jy_{mj}+a_1h^\beta \underset{j=0}{\overset{m}{}}c_jy_{mj}+a_0y_m=u_m$$
(4)
From the approximation (4) we can derive  the following explicit recursive formula for the calculation of the values $`y_m(m=2,3,\mathrm{})`$

$$y_m=\frac{u_ma_2h^\alpha _{j=1}^mb_jy_{mj}a_1h^\beta _{j=1}^mc_jy_{mj}}{a_2h^\alpha b_0+a_1h^\beta c_0+a_0}$$
(5)
with $`y_0=0,y_1=0`$. An analytical solution of fractional-order differential equations is described in detail in , and for closed regulation circuits also in .

## 4. METHOD OF DIFFERENTIATION OF TRANSFER
CHARACTERISTICS

Assume the input of the system (1) with known values $`\alpha `$ and $`\beta `$ is acted upon by arbitrary function $`u(t)`$, and $`y(t)`$ is the output function of the system, while the initial conditions need not be zero in this method. If the input and output functions are known functions of time obtained, for example, by measurement, then for given discrete times $`t_m(m=2,3,\mathrm{})`$ we can determine the values of the derivatives of the output quantity according to the formula (2). For each time instance $`t_m`$ the equation (4) must be satisfied. Therefore after substituting the corresponding values into the equation (4) we obtain a system of three linear equations for unknown coefficients $`a_0,a_1,a_2`$. By solving this system of equations and obtaining the values of coefficients $`a_0,a_1,a_2`$ the system, from which we used experimental data in the calculations, is identified. In order to be able to use more measured values $`y_i,u_i`$, and to make the identification more accurate, we consider for the estimate of the vector of unknown parameters $`\overline{a}`$ as an error criterion the following functional

$$E(\overline{a})=\underset{0}{\overset{T}{}}[a_2y^{(\alpha )}(t)+a_1y^{(\beta )}(t)+a_0y(t)u(t)]^2𝑑tmin$$
(6)
A necessary condition for the minimum to be achieved is

$$\frac{E(\overline{a})}{\overline{a}}=0$$
(7)
Hence, after manipulations we obtain a system of three linear equations for the vector of unknown parameters $`\overline{a}`$

$`a_2{\displaystyle \underset{0}{\overset{T}{}}}(y^{(\alpha )}(t))^2𝑑t+a_1{\displaystyle \underset{0}{\overset{T}{}}}y^{(\alpha )}(t)y^{(\beta )}(t)𝑑t+a_0{\displaystyle \underset{0}{\overset{T}{}}}y^{(\alpha )}(t)y(t)𝑑t={\displaystyle \underset{0}{\overset{T}{}}}y^{(\alpha )}(t)u(t)𝑑t`$
$$a_2\underset{0}{\overset{T}{}}y^{(\alpha )}(t)y^{(\beta )}(t)𝑑t+a_1\underset{0}{\overset{T}{}}(y^{(\beta )}(t))^2𝑑t+a_0\underset{0}{\overset{T}{}}y^{(\beta )}(t)y(t)𝑑t=\underset{0}{\overset{T}{}}y^{(\beta )}(t)u(t)𝑑t$$
(8)
$`a_2{\displaystyle \underset{0}{\overset{T}{}}}y^{(\alpha )}(t)y(t)𝑑t+a_1{\displaystyle \underset{0}{\overset{T}{}}}y^{(\beta )}(t)y(t)𝑑t+a_0{\displaystyle \underset{0}{\overset{T}{}}}(y(t))^2𝑑t={\displaystyle \underset{0}{\overset{T}{}}}y(t)u(t)𝑑t`$
where $`y(t)`$ and $`u(t)`$ are the measured output and input functions of the system, and the derivatives $`y^{(\alpha )}(t)`$ and $`y^{(\beta )}(t)`$ are the derivatives of the output function of the system, obtained by calculation from the formula (2). In the discretization of the system (8) we substitute the continuous functions $`y(t),y^{(\alpha )}(t),y^{(\beta )}(t),`$ $`u(t)`$ in discrete time instants $`t_m(m=0,1,2,3,\mathrm{},M)`$ with their discrete values $`y_m,y_m^{(\alpha )},y_m^{(\beta )},u_m`$ and the integration we substitute with the summation for the time interval $`T`$ for which we have $`M+1`$ measured values from the system. Then we get

$`a_2{\displaystyle \underset{m=0}{\overset{M}{}}}(y_m^{(\alpha )})^2+a_1{\displaystyle \underset{m=0}{\overset{M}{}}}y_m^{(\alpha )}y_m^{(\beta )}+a_0{\displaystyle \underset{m=0}{\overset{M}{}}}y_m^{(\alpha )}y_m={\displaystyle \underset{m=0}{\overset{M}{}}}y_m^{(\alpha )}u_m`$
$$a_2\underset{m=0}{\overset{M}{}}y_m^{(\alpha )}y_m^{(\beta )}+a_1\underset{m=0}{\overset{M}{}}(y_m^{(\beta )})^2+a_0\underset{m=0}{\overset{M}{}}y_m^{(\beta )}y_m=\underset{m=0}{\overset{M}{}}y_m^{(\beta )}u_m$$
(9)
$`a_2{\displaystyle \underset{m=0}{\overset{M}{}}}y_m^{(\alpha )}y_m+a_1{\displaystyle \underset{m=0}{\overset{M}{}}}y_m^{(\beta )}y_m+a_0{\displaystyle \underset{m=0}{\overset{M}{}}}(y_m)^2={\displaystyle \underset{m=0}{\overset{M}{}}}y_mu_m`$
From the system of equations (9) we can calculate the unknown parameters $`\overline{a}`$, and the task of parameter identification of the system is thus solved provided we know the values $`\alpha `$ and $`\beta `$.

## 5. METHOD OF IDENTIFICATION OF FRACTIONAL-ORDER
DYNAMICAL SYSTEMS

In general in the sytems of the type (1) it is necessary to identify the orders of derivatives $`\alpha `$ and $`\beta `$ which are non-integer and real. If, under this assumption, we used the given method of differentiation of transfer characteristics for the system (1), we would obtain a system of five nonlinear equations. It is difficult to solve such systems of equations, however, therefore we chose a combination of the method of differentiation of transfer characteristics for the calculation of the coefficients $`a_0,a_1,a_2`$ and the method of simple search for the determination of the coefficients $`\alpha `$ and $`\beta `$. While the first method is analytical, in the second the coefficients are calculated iteratively.

In the identification of the parameters of the model we proceed as follows:

For the parameters $`\alpha `$ and $`\beta `$ we determine the intervals $`<\alpha _{min},\alpha _{max}>`$ and $`<\beta _{min},\beta _{max}>`$, in which we look for the parameters. We divide each of these intervals into $`n_\alpha `$ and $`n_\beta `$ subintervals

$$n_\alpha =\frac{2(\alpha _{max}\alpha _{min})}{\epsilon }$$
(10)
$$n_\beta =\frac{2(\beta _{max}\beta _{min})}{\epsilon }$$
(11)
where $`\epsilon `$ denotes the fineness of division. We choose the starting values $`\alpha `$ and $`\beta `$ from the chosen intervals.

We use the method of differentiation of the transfer characteristic to calculate for given $`\alpha `$ and $`\beta `$ the corresponding values of the parameters $`a_0,a_1,a_2`$.

We calculate the criterion of the approximation of experimental output values of the real system $`ye_i`$ and the outputs of the model $`y_i`$

$$Q=\frac{1}{M+1}\underset{i=0}{\overset{M}{}}[ye_iy_i]^2$$
(12)
which serves for the choice of optimum values of the coefficients $`\alpha `$ and $`\beta `$.

If the value of the actually computed approximation criterion is less than the value of the criterion from the preceding iteration step, we choose for the optimum subintervals those containing last values of the coefficients $`\alpha `$ and $`\beta `$.

After the calculation of the approximation criterion for all subintervals, we determine from the optimum subintervals the intervals $`<\alpha _{min},\alpha _{max}>`$ and $`<\beta _{min},\beta _{max}>`$ and for these subintervals the new subintervals. We repeat this procedure until the required accuracy of the calculation of parameters $`\alpha `$ and $`\beta `$ is achieved, where we have always calculated the values of the parameters $`a_0,a_1,a_2`$ by using the method of differentiation of transfer characteristics.

## 6. VERIFICATION OF THE METHOD

The verification of the correctness of the described method of parameter identification will be done by using it in systems with known parameters.

Consider an integer-order system  with coefficients $`a_2=1,a_1=3,a_0=2,\alpha =2,\beta =1.`$ The transfer characteristic of this system has aperiodic behavior and is shown in Fig.1. For the chosen intervals $`\alpha <2;2>`$ and $`\beta <1;1>`$ the precise values of the coefficients $`a_2=1,a_1=3,a_0=2`$ were calculated with the method of differentiation of transfer characteristics. For the intervals $`\alpha <1,5;2,55>`$ and $`\beta <0,7;1,33>`$ the following values of the coefficients were calculated with both methods $`a_2=1,0001,`$ $`a_1=2,99987,a_0=1,99998,\alpha =1,99993,\beta =0,99996`$ , which values are very close to the parameters of the identified system.

Consider a fractional-order system with coefficients  $`a_2=0,8,a_1=0,5,a_0=1,\alpha =2,2,\beta =0,9`$. The transfer characteristic of this system computed according to the relation (5), has periodic behavior and is shown in Fig.2. For the chosen intervals $`\alpha <2,2;2,2>`$ and $`\beta <0,9;0,9>`$ the precise values of the coefficients $`a_2=0,8,a_1=0,5,a_0=1`$ were computed with the method of differentiation of transfer characteristics. For the intervals $`\alpha <1,5;2,55>`$ and $`\beta <0,7;1,33>`$ the following values were computed by using both methods $`a_2=0,80005,`$ $`a_1=0,49996,a_0=0,99998,\alpha =2,19996,\beta =0,89989`$ , which values are very close to the parameters of the identified system.

Let us consider again the previous fractional-order system with coefficients $`a_2=0,8,`$ $`a_1=0,5,a_0=1,\alpha =2,2,\beta =0,9,`$ whose behavior is given in Fig.2. Let us identify it as a second-order system, which was the classical procedure up until now. If we want to identify the system as a second-order system, we have to specify the same upper and lower bounds $`\alpha <2;2>`$ and $`\beta <1;1>.`$ With the method of differentiation of transfer characteristics the following coefficient values were calculated : $`a_2=0,76639,a_1=0,23184,a_0=1.`$ The system with these parameters is actually a classical integer-order model of a fractional-order system; their graphical comparison is shown in Fig.3. Although optically both characteristics are rather close to each other, the use of such model in the synthesis of the parameters of the regulator is associated with unfavorable consequences for the dynamics of the regulation circuit as shown e.g. in .

Let us now identify the parameters of an experimental heating furnace whose transfer characteristic is shown in Fig. 4. Assuming it is a second order system with $`\alpha =2,\beta =1,`$ we obtain parameters $`a_0=1,928,a_1=4892,733,a_2=73043,36`$ and the value $`1,02.10^3`$ of the criterion (12). For the chosen intervals $`\alpha <1,1;2,55>`$ and $`\beta <0,33;1,3>`$ , the parameters $`a_2=14994,3,`$ $`a_1=6009,52,a_0=1,69,\alpha =1,31,\beta =0,97`$ were identified for this system, and the criterion (12) had the value $`2,7.10^4`$ (see Fig. 4). As we can see ($`a_2`$) the previous two approximations are unsuitable. Assuming a two-member differential equation, the parameters take on the values $`a_1=788,35,a_0=1,39,\beta =0,73`$ and the value of the criterion (12) is $`6,3.10^4`$.

## 7. CONCLUSION

The above results of parameter identification for fractional- or integer-order dynamical systems confirm that in known systems the agreement of calculated and actual parameters was very good. In the case of real objects identification it is necessary to pay attention to the accuracy of the output quantity measurement in the identified system, and to the time step of the measurement from the point of view of algorithms.

When identifying the parameters of fractionarl-oder dynamical systems in which real orders of derivatives $`\alpha `$ and $`\beta `$ are identified, it is suitable, considering the existence of local minima of the functional (6) , to repeat the calculations for different intervals $`<\alpha _{min},\alpha _{max}>`$ and $`<\beta _{min},\beta _{max}>`$, in which we look for the parameters, and for different fineness of division $`\epsilon `$ in the formulas (10) and (11).

## 8. REFERENCES

 K. B. Oldham and J. Spanier: The Fractional Calculus. Academic Press, New York, 1974.

 I. Podlubny : The Laplace Transform Method for Linear Differential Equations of the Fractional Order, UEF - 02 - 94, The Academy of Sciences Institute of Experimental Physics, 1994, Kosice.

 M. Axtell and M. E. Bise : Fractional Calculus Applications in Control Systems, in: Proc. of the IEEE 1990 Nat. Aerospace and Electronics Conf., New York, 1990, pp. 563 - 566.

 Ľ. Dorčák : Numerical Models for Simulation of the Fractional-Order Control Systems, in: UEF - 04 - 94, The Academy of Sciences Institute of Experimental Physics, 1994, Kosice, p. 12.

 Ľ. Sýkorová, I. Koštial, Ľ. Dorčák : A comparison of Integer- and Fractional- Order PID Regulators, in: Proceeding of the 2nd scientific-technical conference with international participation PROCESS CONTROL, Horná Bečva, Czech Republic, June 3-6, 1996, pp. 347-351.

 I. Podlubny. Numerical methods of the fractional calculus, part I, II. Trans. of the TU Kosice, vol. 4, no. 3-4, pp. 200-208.