function [phasevel,atten,groupvel]=cal_CQU(freq)

% function designed to interpolate phase velocity and Q
% in PREM model by using Jean's formular and linear interpolation
% Usage: [a,b]=cal_CQ(freq);
% a: phase velocoty vector
% b: attenuation vector
% Using prem228 model (from Yuancheng)

% group velocity and attenuation w.r.t period in PREM model.
% angular order, period, group velocity and attenuation factor(Q)
PREM = [ ...
2 4119.604 5.989203 310.4874;
3 2462.702 7.139105 320.8211;
4 1677.931 7.950159 326.2403;
5 1250.607 8.274653 329.8324;
6 992.8075 8.287525 331.4599;
7 825.1759 8.031121 330.5515;
8 711.0776 7.467957 326.5022;
9 631.6581 6.650135 319.7901;
10 574.9651 5.867804 312.3260;
11 532.2573 5.336467 305.0800;
12 498.0367 5.013240 297.7057;
13 469.3875 4.799463 289.8472;
14 444.7692 4.636135 281.4356;
15 423.2746 4.496884 272.5742;
16 404.2994 4.371479 263.4450;
17 387.4018 4.256602 254.2537;
18 372.2400 4.151586 245.1935;
19 358.5410 4.056571 236.4244;
20 346.0831 3.971742 228.0647;
21 334.6846 3.897062 220.1904;
22 324.1957 3.832229 212.8412;
23 314.4925 3.776710 206.0278;
24 305.4721 3.729807 199.7406;
25 297.0491 3.690728 193.9558;
26 289.1521 3.658644 188.6416;
27 281.7211 3.632729 183.7621;
28 274.7060 3.612192 179.2800;
29 268.0641 3.596299 175.1586;
30 261.7592 3.584379 171.3631;
31 255.7606 3.575829 167.8611;
32 250.0416 3.570118 164.6229;
33 244.5793 3.566779 161.6214;
34 239.3537 3.565407 158.8325;
35 234.3473 3.565652 156.2344;
36 229.5447 3.567216 153.8079;
37 224.9321 3.569845 151.5356;
38 220.4973 3.573322 149.4025;
39 216.2295 3.577466 147.3950;
40 212.1186 3.582125 145.5011;
41 208.1557 3.587170 143.7102;
42 204.3327 3.592498 142.0130;
43 200.6419 3.598019 140.4011;
44 197.0766 3.603662 138.8671;
45 193.6304 3.609369 137.4046;
46 190.2973 3.615092 136.0078;
47 187.0720 3.620792 134.6714;
48 183.9493 3.626440 133.3909;
49 180.9244 3.632011 132.1623;
50 177.9930 3.637488 130.9818;
51 175.1508 3.642855 129.8463;
52 172.3939 3.648104 128.7528;
53 169.7187 3.653227 127.6987;
54 167.1217 3.658219 126.6818;
55 164.5995 3.663078 125.6997;
56 162.1491 3.667804 124.7508;
57 159.7676 3.672396 123.8332;
58 157.4521 3.676856 122.9454;
59 155.2001 3.681187 122.0860;
60 153.0091 3.685391 121.2537;
61 150.8766 3.689472 120.4473;
62 148.8006 3.693434 119.6658;
63 146.7787 3.697281 118.9080;
64 144.8090 3.701016 118.1732;
65 142.8895 3.704645 117.4604;
66 141.0185 3.708170 116.7689;
67 139.1940 3.711598 116.0979;
68 137.4146 3.714931 115.4467;
69 135.6785 3.718174 114.8148;
70 133.9843 3.721331 114.2014;
71 132.3305 3.724406 113.6059;
72 130.7156 3.727402 113.0280;
73 129.1385 3.730323 112.4669;
74 127.5977 3.733173 111.9223;
75 126.0921 3.735954 111.3936;
76 124.6205 3.738670 110.8805;
77 123.1819 3.741323 110.3824;
78 121.7751 3.743917 109.8989;
79 120.3990 3.746454 109.4297;
80 119.0529 3.748936 108.9744;
81 117.7356 3.751366 108.5326;
82 116.4463 3.753746 108.1038;
83 115.1841 3.756078 107.6879;
84 113.9482 3.758364 107.2844;
85 112.7379 3.760606 106.8930;
86 111.5522 3.762805 106.5134;
87 110.3906 3.764963 106.1452;
88 109.2522 3.767081 105.7882;
89 108.1364 3.769161 105.4422;
90 107.0426 3.771204 105.1067;
91 105.9701 3.773210 104.7815;
92 104.9183 3.775182 104.4664;
93 103.8866 3.777120 104.1610;
94 102.8745 3.779025 103.8652;
95 101.8815 3.780897 103.5787;
96 100.9069 3.782739 103.3012;
97 99.95030 3.784549 103.0325;
98 99.01122 3.786329 102.7723;
99 98.08919 3.788080 102.5205;
100 97.18374 3.789801 102.2768;
101 96.29445 3.791494 102.0410;
102 95.42088 3.793159 101.8129;
103 94.56264 3.794796 101.5923;
104 93.71932 3.796406 101.3790;
105 92.89055 3.797988 101.1727;
106 92.07597 3.799544 100.9734;
107 91.27521 3.801074 100.7808;
108 90.48794 3.802577 100.5947;
109 89.71382 3.804054 100.4150;
110 88.95252 3.805505 100.2414;
111 88.20375 3.806931 100.0739;
112 87.46719 3.808331 99.91229;
113 86.74256 3.809705 99.75636;
114 86.02957 3.811055 99.60597;
115 85.32795 3.812379 99.46098;
116 84.63744 3.813679 99.32124;
117 83.95777 3.814953 99.18662;
118 83.28870 3.816203 99.05696;
119 82.62999 3.817428 98.93214;
120 81.98140 3.818629 98.81203;
121 81.34270 3.819805 98.69649;
122 80.71368 3.820957 98.58540;
123 80.09412 3.822085 98.47865;
124 79.48382 3.823188 98.37610;
125 78.88256 3.824267 98.27765;
126 78.29016 3.825323 98.18318;
127 77.70642 3.826354 98.09258;
128 77.13116 3.827362 98.00575;
129 76.56420 3.828346 97.92259;
130 76.00536 3.829306 97.84299;
131 75.45448 3.830243 97.76685;
132 74.91139 3.831156 97.69409;
133 74.37592 3.832046 97.62460;
134 73.84793 3.832913 97.55829;
135 73.32726 3.833757 97.49509;
136 72.81376 3.834578 97.43489;
137 72.30729 3.835376 97.37763;
138 71.80770 3.836152 97.32321;
139 71.31487 3.836905 97.27157;
140 70.82865 3.837635 97.22261;
141 70.34893 3.838344 97.17628;
142 69.87556 3.839030 97.13249;
143 69.40844 3.839694 97.09117;
144 68.94743 3.840337 97.05227;
145 68.49243 3.840958 97.01570;
146 68.04331 3.841558 96.98141;
147 67.59998 3.842136 96.94932;
148 67.16231 3.842693 96.91939;
149 66.73020 3.843230 96.89155;
150 66.30356 3.843746 96.86574;
151 65.88228 3.844241 96.84191;
152 65.46625 3.844716 96.81999;
153 65.05540 3.845171 96.79995;
154 64.64962 3.845606 96.78172;
155 64.24881 3.846021 96.76525;
156 63.85291 3.846416 96.75050;
157 63.46180 3.846793 96.73742;
158 63.07542 3.847150 96.72596;
159 62.69367 3.847488 96.71607;
160 62.31649 3.847808 96.70772;
161 61.94377 3.848109 96.70085;
162 61.57546 3.848392 96.69543;
163 61.21148 3.848658 96.69142;
164 60.85174 3.848905 96.68878;
165 60.49618 3.849134 96.68747;
166 60.14473 3.849347 96.68745;
167 59.79731 3.849542 96.68869;
168 59.45387 3.849720 96.69114;
169 59.11433 3.849882 96.69479;
170 58.77863 3.850027 96.69958;
171 58.44671 3.850156 96.70550;
172 58.11850 3.850269 96.71251;
173 57.79395 3.850366 96.72058;
174 57.47299 3.850448 96.72968;
175 57.15557 3.850514 96.73978;
176 56.84163 3.850565 96.75086;
177 56.53111 3.850602 96.76288;
178 56.22396 3.850623 96.77581;
179 55.92013 3.850631 96.78965;
180 55.61957 3.850624 96.80435;
181 55.32221 3.850603 96.81989;
182 55.02803 3.850569 96.83626;
183 54.73695 3.850520 96.85342;
184 54.44894 3.850459 96.87136;
185 54.16395 3.850385 96.89005;
186 53.88194 3.850297 96.90947;
187 53.60285 3.850197 96.92961;
188 53.32664 3.850085 96.95043;
189 53.05328 3.849960 96.97193;
190 52.78271 3.849824 96.99409;
191 52.51490 3.849675 97.01688;
192 52.24980 3.849515 97.04028;
193 51.98737 3.849344 97.06429;
194 51.72759 3.849161 97.08888;
195 51.47040 3.848967 97.11404;
196 51.21576 3.848763 97.13975;
197 50.96365 3.848548 97.16600;
198 50.71403 3.848323 97.19277;
199 50.46685 3.848087 97.22005;
200 50.22209 3.847841 97.24782];

Erad = 6371; % Radius of the Earth in km
Pvel=2*pi./PREM(:,2).*Erad./(PREM(:,1)+0.5); % by Jean's formular
Q=PREM(:,4);
Tuc=PREM(:,3);

% find out an appropriate interval.

tfreq=1./PREM(:,2);
minfreq=min(tfreq);
maxfreq=max(tfreq);
nmod = size(PREM,1);

nn=length(freq);
phasevel=zeros(nn,1);
atten=zeros(nn,1);
groupvel=zeros(nn,1);

for j=1:nn

if (freq(j) < minfreq) | ( freq(j) > maxfreq)
   if (freq(j) > maxfreq)
%      phasevel(j)=Pvel(nmod);

      Ftmp2=(log10(freq(j))-log10(tfreq(nmod)))/ ...
            (log10(tfreq(nmod))-log10(tfreq(nmod-1)));
      atten(j)=Ftmp2*(Q(nmod)-Q(nmod-1))+Q(nmod);
      phasevel(j)=Ftmp2*(Pvel(nmod)-Pvel(nmod-1))+Pvel(nmod);
      groupvel(j)=Ftmp2*(Tuc(nmod)-Tuc(nmod-1))+Tuc(nmod);
   end
   if (freq(j) < minfreq)
      if (freq(j) == 0)
         freq(j)=freq(j+1);
      end
      Ftmp1=(log10(tfreq(1))-log10(freq(j)))/ ...
            (log10(tfreq(2))-log10(freq(j))); 
      phasevel(j)=(Ftmp1*Pvel(2)-Pvel(1))/(Ftmp1-1);
      atten(j)=(Ftmp1*Q(2)-Q(1))/(Ftmp1-1);    
      groupvel(j)=(Ftmp1*Tuc(2)-Tuc(1))/(Ftmp1-1);
        
   end
else
   for i=1:nmod-1
      if (freq(j) >= tfreq(i)) & (freq(j) <= tfreq(i+1))
          Ptmp=Pvel(i+1)-Pvel(i);
          Gtmp=Tuc(i+1)-Tuc(i);
          Qtmp=Q(i+1)-Q(i);
          Ftmp=(log10(freq(j))-log10(tfreq(i)))/ ...
               (log10(tfreq(i+1))-log10(tfreq(i)));  
          phasevel(j) = Ptmp*Ftmp+Pvel(i);
          atten(j) = Qtmp*Ftmp+Q(i);
          groupvel(j) = Gtmp*Ftmp+Tuc(i);
      end
   end
end    

end

