source/unres/src_MD-M/intlocal.f

   1       subroutine integral(gamma1,gamma2,gamma3,gamma4,ity1,ity2,a1,a2,
   2      &  si1,si2,si3,si4,transp,q)
   3       implicit none
   4       integer ity1,ity2
   5       integer ilam1,ilam2,ilam3,ilam4,iincr
   6       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
   7      &  lambda2,lambda3,lambda4
   8       logical transp
   9       double precision elocal,ele
  10       double precision delta,delta2,sum,ene,sumene,boltz
  11       double precision q,a1(2,2),a2(2,2),si1,si2,si3,si4
  12       double precision conv /.01745329252d0/,pi /3.141592654d0/
  13
  14       iincr=20
  15       delta=iincr*conv
  16       delta2=0.5d0*delta
  17 cd      print *,'iincr',iincr,' delta',delta
  18 cd      write(2,*) gamma1,gamma2,ity1,ity2,a1,a2,si1,si2,si3,si4,transp
  19
  20 cd      do ilam1=-180,180,5
  21 cd        do ilam2=-180,180,5
  22 cd          lambda1=ilam1*conv+delta2
  23 cd          lambda2=ilam2*conv+delta2
  24 cd          write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
  25 cd     &    ele(lambda1,lambda2,a1,1.0d0,1.d00)
  26 cd        enddo
  27 cd      enddo
  28 cd      stop
  29
  30       sum=0.0d0
  31       sumene=0.0d0
  32       do ilam1=-180,179,iincr
  33         do ilam2=-180,179,iincr
  34           do ilam3=-180,179,iincr
  35             do ilam4=-180,179,iincr
  36               lambda1=ilam1*conv+delta2
  37               lambda2=ilam2*conv+delta2
  38               lambda3=ilam3*conv+delta2
  39               lambda4=ilam4*conv+delta2
  40 cd              write (2,*) ilam1,ilam2,ilam3,ilam4
  41 cd              write (2,*) lambda1,lambda2,lambda3,lambda4
  42               ene=
  43      &         -elocal(ity1,lambda1,lambda2,.false.)*
  44      &          elocal(ity2,lambda3,lambda4,transp)*
  45      &          ele(si1*lambda1+gamma1,si3*lambda3+gamma3,a1)*
  46      &          ele(si2*lambda2+gamma2,si4*lambda4+gamma4,a2)
  47 cd              write (2,*) elocal(ity1,lambda1,gamma1-pi-lambda2),
  48 cd     &        elocal(ity2,lambda3,gamma2-pi-lambda4),
  49 cd     &        ele(lambda1,lambda2,a1,si1,si3),
  50 cd     &        ele(lambda3,lambda4,a2,si2,si4)
  51               sum=sum+ene
  52             enddo
  53           enddo
  54         enddo
  55       enddo
  56       q=sum/(2*pi)**4*delta**4
  57       write (2,* )'sum',sum,' q',q
  58       return
  59       end
  60 c---------------------------------------------------------------------------
  61       subroutine integral3(gamma1,gamma2,ity1,ity2,ity3,ity4,
  62      &  a1,koniec,q1,q2,q3,q4)
  63       implicit none
  64       integer ity1,ity2,ity3,ity4
  65       integer ilam1,ilam2,ilam3,ilam4,iincr
  66       double precision gamma1,gamma2,gamma3,gamma4,beta,lambda1,
  67      &  lambda2,lambda3,lambda4
  68       logical koniec
  69       double precision elocal,ele
  70       double precision delta,delta2,sum1,sum2,sum3,sum4,
  71      &  ene1,ene2,ene3,ene4,boltz
  72       double precision q1,q2,q3,q4,a1(2,2),a2(2,2)
  73       double precision conv /.01745329252d0/,pi /3.141592654d0/
  74
  75       iincr=60
  76       delta=iincr*conv
  77       delta2=0.5d0*delta
  78 cd      print *,'iincr',iincr,' delta',delta
  79       write(2,*) gamma1,gamma2,ity1,ity2,ity3,ity4,a1,koniec
  80
  81 cd      do ilam1=-180,180,5
  82 cd        do ilam2=-180,180,5
  83 cd          lambda1=ilam1*conv+delta2
  84 cd          lambda2=ilam2*conv+delta2
  85 cd          write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
  86 cd     &    ele(lambda1,lambda2,a1,1.0d0,1.d00)
  87 cd        enddo
  88 cd      enddo
  89 cd      stop
  90
  91       sum1=0.0d0
  92       sum2=0.0d0
  93       sum3=0.0d0
  94       sum4=0.0d0
  95       do ilam1=-180,179,iincr
  96         do ilam2=-180,179,iincr
  97           do ilam3=-180,179,iincr
  98             do ilam4=-180,179,iincr
  99               lambda1=ilam1*conv+delta2
 100               lambda2=ilam2*conv+delta2
 101               lambda3=ilam3*conv+delta2
 102               lambda4=ilam4*conv+delta2
 103 cd              write (2,*) ilam1,ilam2,ilam3,ilam4
 104 cd              write (2,*) lambda1,lambda2,lambda3,lambda4
 105               if (.not.koniec) then
 106               ene1=
 107      &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
 108      &          elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
 109      &          ele(lambda2,lambda4,a1)
 110               else
 111               ene1=
 112      &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
 113      &          elocal(ity3,lambda3,lambda4,.false.)*
 114      &          ele(lambda2,-lambda4,a1)
 115               endif
 116               ene2=
 117      &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
 118      &          elocal(ity4,lambda3,lambda4,.false.)*
 119      &          ele(lambda2,lambda3,a1)
 120               if (.not.koniec) then
 121               ene3=
 122      &          elocal(ity2,lambda1,lambda2,.false.)*
 123      &          elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
 124      &          ele(lambda1,lambda4,a1)
 125               else
 126               ene3=
 127      &          elocal(ity2,lambda1,lambda2,.false.)*
 128      &          elocal(ity3,lambda3,lambda4,.false.)*
 129      &          ele(lambda1,-lambda4,a1)
 130               endif
 131               ene4=
 132      &          elocal(ity2,lambda1,lambda2,.false.)*
 133      &          elocal(ity4,lambda3,lambda4,.false.)*
 134      &          ele(lambda1,lambda3,a1)
 135               sum1=sum1+ene1
 136               sum2=sum2+ene2
 137               sum3=sum3+ene3
 138               sum4=sum4+ene4
 139             enddo
 140           enddo
 141         enddo
 142       enddo
 143       q1=sum1/(2*pi)**4*delta**4
 144       q2=sum2/(2*pi)**4*delta**4
 145       q3=sum3/(2*pi)**4*delta**4
 146       q4=sum4/(2*pi)**4*delta**4
 147       write (2,* )'sum',sum1,sum2,sum3,sum4,' q',q1,q2,q3,q4
 148       return
 149       end
 150 c-------------------------------------------------------------------------
 151       subroutine integral5(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
 152      &  ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4)
 153       implicit none
 154       integer ity1,ity2,ity3,ity4,ity5,ity6
 155       integer ilam1,ilam2,ilam3,ilam4,ilam5,iincr
 156       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
 157      &  lambda2,lambda3,lambda4,lambda5
 158       logical transp
 159       double precision elocal,ele
 160       double precision eloc1,eloc2,eloc3,eloc4,eloc5,eloc6,ele1,ele2
 161       double precision delta,delta2,sum,ene,sumene,pom
 162       double precision ene1,ene2,ene3,ene4,sum1,sum2,sum3,sum4,
 163      &  a1(2,2),a2(2,2)
 164       integer si1,si2,si3,si4
 165       double precision conv /.01745329252d0/,pi /3.141592654d0/
 166
 167       iincr=60
 168       delta=iincr*conv
 169       delta2=0.5d0*delta
 170 cd      print *,'iincr',iincr,' delta',delta
 171 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
 172 cd     &  ' gamma3=',gamma3,' gamma4=',gamma4
 173 cd      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
 174 cd      write(2,*) 'a1=',a1
 175 cd      write(2,*) 'a2=',a2
 176 cd      write(2,*) si1,si2,si3,si4,transp
 177
 178       sum1=0.0d0
 179       sum2=0.0d0
 180       sum3=0.0d0
 181       sum4=0.0d0
 182       do ilam1=-180,179,iincr
 183         do ilam2=-180,179,iincr
 184           do ilam3=-180,179,iincr
 185             do ilam4=-180,179,iincr
 186               do ilam5=-180,179,iincr
 187                 lambda1=ilam1*conv+delta2
 188                 lambda2=ilam2*conv+delta2
 189                 lambda3=ilam3*conv+delta2
 190                 lambda4=ilam4*conv+delta2
 191                 lambda5=ilam5*conv+delta2
 192                 if (transp) then
 193                   ele1=ele(lambda1,si4*lambda4,a1)
 194                   ele2=ele(lambda2,lambda3,a2)
 195                 else
 196                   ele1=ele(lambda1,lambda3,a1)
 197                   ele2=ele(lambda2,si4*lambda4,a2)
 198                 endif
 199                 eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
 200                 eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
 201                 pom=ele1*ele2*eloc2*eloc5
 202                 if (si1.gt.0) then
 203                   eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
 204                   sum1=sum1+pom*eloc1
 205                 endif
 206                 eloc3=elocal(ity3,lambda2,lambda5,.false.)
 207                 sum2=sum2+pom*eloc3
 208                 eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
 209                 sum3=sum3+pom*eloc4
 210                 if (si4.gt.0) then
 211                   eloc6=elocal(ity6,lambda4,lambda5,.false.)
 212                   sum4=sum4+pom*eloc6
 213                 endif
 214               enddo
 215             enddo
 216           enddo
 217         enddo
 218       enddo
 219       pom=1.0d0/(2*pi)**5*delta**5
 220       ene1=sum1*pom
 221       ene2=sum2*pom
 222       ene3=sum3*pom
 223       ene4=sum4*pom
 224 c      write (2,* )'sum',sum1,sum2,sum3,sum4,' q',ene1,ene2,ene3,ene4
 225       return
 226       end
 227 c-------------------------------------------------------------------------
 228       subroutine integral_turn6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,
 229      &  ity3,ity4,ity5,ity6,a1,a2,ene_turn6)
 230       implicit none
 231       integer ity1,ity2,ity3,ity4,ity5,ity6
 232       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
 233       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
 234      &  lambda2,lambda3,lambda4,lambda5,lambda6
 235       logical transp
 236       double precision elocal,ele
 237       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
 238      &  eloc61,ele1,ele2
 239       double precision delta,delta2,sum,ene,sumene,pom,ene5
 240       double precision ene_turn6,sum5,a1(2,2),a2(2,2)
 241       double precision conv /.01745329252d0/,pi /3.141592654d0/
 242
 243       iincr=60
 244       delta=iincr*conv
 245       delta2=0.5d0*delta
 246 cd      print *,'iincr',iincr,' delta',delta
 247       write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
 248      &  ' gamma3=',gamma3,' gamma4=',gamma4
 249       write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
 250       write(2,*) 'a1=',a1
 251       write(2,*) 'a2=',a2
 252
 253       sum5=0.0d0
 254       do ilam1=-180,179,iincr
 255         do ilam2=-180,179,iincr
 256           do ilam3=-180,179,iincr
 257             do ilam4=-180,179,iincr
 258               do ilam5=-180,179,iincr
 259                 lambda1=ilam1*conv+delta2
 260                 lambda2=ilam2*conv+delta2
 261                 lambda3=ilam3*conv+delta2
 262                 lambda4=ilam4*conv+delta2
 263                 lambda5=ilam5*conv+delta2
 264                 ele1=ele(lambda1,-lambda4,a1)
 265                 ele2=ele(lambda2,lambda3,a2)
 266                 eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
 267                 eloc5=elocal(ity5,lambda3,lambda4,.false.)
 268                 pom=ele1*ele2*eloc2*eloc5
 269                 eloc3=elocal(ity3,lambda2,gamma3-pi-lambda5,.false.)
 270                 eloc4=elocal(ity4,lambda5,gamma4-pi-lambda3,.false.)
 271                 sum5=sum5+pom*eloc3*eloc4
 272               enddo
 273             enddo
 274           enddo
 275         enddo
 276       enddo
 277       pom=-1.0d0/(2*pi)**5*delta**5
 278       ene_turn6=sum5*pom
 279 c      print *,'sum6',sum6,' ene6',ene6
 280       return
 281       end
 282 c-------------------------------------------------------------------------
 283       subroutine integral6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
 284      &  ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4,
 285      &  ene5,ene6)
 286       implicit none
 287       integer ity1,ity2,ity3,ity4,ity5,ity6
 288       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
 289       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
 290      &  lambda2,lambda3,lambda4,lambda5,lambda6
 291       logical transp
 292       double precision elocal,ele
 293       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
 294      &  eloc61,ele1,ele2
 295       double precision delta,delta2,sum,ene,sumene,pom
 296       double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
 297      &  sum4,sum5,sum6,a1(2,2),a2(2,2)
 298       integer si1,si2,si3,si4
 299       double precision conv /.01745329252d0/,pi /3.141592654d0/
 300
 301       iincr=60
 302       delta=iincr*conv
 303       delta2=0.5d0*delta
 304 cd      print *,'iincr',iincr,' delta',delta
 305 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
 306 cd     &  ' gamma3=',gamma3,' gamma4=',gamma4
 307 cd      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
 308 cd      write(2,*) 'a1=',a1
 309 cd      write(2,*) 'a2=',a2
 310 cd      write(2,*) si1,si2,si3,si4,transp
 311
 312       sum1=0.0d0
 313       sum2=0.0d0
 314       sum3=0.0d0
 315       sum4=0.0d0
 316       sum5=0.0d0
 317       sum6=0.0d0
 318       eloc1=0.0d0
 319       eloc6=0.0d0
 320       eloc61=0.0d0
 321       do ilam1=-180,179,iincr
 322         do ilam2=-180,179,iincr
 323           do ilam3=-180,179,iincr
 324             do ilam4=-180,179,iincr
 325               do ilam5=-180,179,iincr
 326                 do ilam6=-180,179,iincr
 327                 lambda1=ilam1*conv+delta2
 328                 lambda2=ilam2*conv+delta2
 329                 lambda3=ilam3*conv+delta2
 330                 lambda4=ilam4*conv+delta2
 331                 lambda5=ilam5*conv+delta2
 332                 lambda6=ilam6*conv+delta2
 333                 if (transp) then
 334                   ele1=ele(lambda1,si4*lambda4,a1)
 335                   ele2=ele(lambda2,lambda3,a2)
 336                 else
 337                   ele1=ele(lambda1,lambda3,a1)
 338                   ele2=ele(lambda2,si4*lambda4,a2)
 339                 endif
 340                 eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
 341                 eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
 342                 pom=ele1*ele2*eloc2*eloc5
 343                 if (si1.gt.0) then
 344                   eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
 345                 endif
 346                 eloc3=elocal(ity3,lambda2,lambda6,.false.)
 347                 sum1=sum1+pom*eloc1*eloc3
 348                 eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
 349                 if (si4.gt.0) then
 350                   eloc6=elocal(ity6,lambda4,lambda6,.false.)
 351                   eloc61=elocal(ity6,lambda4,lambda5,.false.)
 352                 endif
 353                 sum2=sum2+pom*eloc4*eloc6
 354                 eloc41=elocal(ity4,lambda6,gamma3-pi-lambda3,.false.)
 355                 sum3=sum3+pom*eloc1*eloc41
 356                 sum4=sum4+pom*eloc1*eloc6
 357                 sum5=sum5+pom*eloc3*eloc4
 358                 sum6=sum6+pom*eloc3*eloc61
 359                 enddo
 360               enddo
 361             enddo
 362           enddo
 363         enddo
 364       enddo
 365       pom=-1.0d0/(2*pi)**6*delta**6
 366       ene1=sum1*pom
 367       ene2=sum2*pom
 368       ene3=sum3*pom
 369       ene4=sum4*pom
 370       ene5=sum5*pom
 371       ene6=sum6*pom
 372 c      print *,'sum6',sum6,' ene6',ene6
 373       return
 374       end
 375 c-------------------------------------------------------------------------
 376       subroutine integral3a(gamma1,gamma2,ity1,ity2,a1,si1,ene1)
 377       implicit none
 378       integer ity1,ity2,ity3,ity4,ity5,ity6
 379       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
 380       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
 381      &  lambda2,lambda3,lambda4,lambda5,lambda6
 382       logical transp
 383       double precision elocal,ele
 384       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
 385      &  eloc61,ele1,ele2
 386       double precision delta,delta2,sum,ene,sumene,pom
 387       double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
 388      &  sum4,sum5,sum6,a1(2,2),a2(2,2)
 389       integer si1,si2,si3,si4
 390       double precision conv /.01745329252d0/,pi /3.141592654d0/
 391
 392       iincr=60
 393       delta=iincr*conv
 394       delta2=0.5d0*delta
 395 cd      print *,'iincr',iincr,' delta',delta
 396 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2
 397 cd      write(2,*) ity1,ity2
 398 cd      write(2,*) 'a1=',a1
 399 cd      write(2,*) si1,
 400
 401       sum1=0.0d0
 402       eloc1=0.0d0
 403       do ilam1=-180,179,iincr
 404         do ilam2=-180,179,iincr
 405           do ilam3=-180,179,iincr
 406             lambda1=ilam1*conv+delta2
 407             lambda2=ilam2*conv+delta2
 408             lambda3=ilam3*conv+delta2
 409             ele1=ele(lambda1,si1*lambda3,a1)
 410             eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
 411             if (si1.gt.0) then
 412               eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
 413             else
 414               eloc2=elocal(ity2,lambda2,lambda3,.false.)
 415             endif
 416             sum1=sum1+ele1*eloc1*eloc2
 417           enddo
 418         enddo
 419       enddo
 420       pom=1.0d0/(2*pi)**3*delta**3
 421       ene1=sum1*pom
 422       return
 423       end
 424 c-------------------------------------------------------------------------
 425       subroutine integral4a(gamma1,gamma2,gamma3,ity1,ity2,ity3,a1,si1,
 426      &  ene1)
 427       implicit none
 428       integer ity1,ity2,ity3,ity4,ity5,ity6
 429       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
 430       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
 431      &  lambda2,lambda3,lambda4,lambda5,lambda6
 432       logical transp
 433       double precision elocal,ele
 434       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
 435      &  eloc61,ele1,ele2
 436       double precision delta,delta2,sum,ene,sumene,pom
 437       double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
 438      &  sum4,sum5,sum6,a1(2,2),a2(2,2)
 439       integer si1,si2,si3,si4
 440       double precision conv /.01745329252d0/,pi /3.141592654d0/
 441
 442       iincr=60
 443       delta=iincr*conv
 444       delta2=0.5d0*delta
 445 cd      print *,'iincr',iincr,' delta',delta
 446 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
 447 cd     &  ' gamma3=',gamma3
 448 cd      write(2,*) ity1,ity2,ity3
 449 cd      write(2,*) 'a1=',a1
 450 cd      write(2,*) 'si1=',si1
 451       sum1=0.0d0
 452       do ilam1=-180,179,iincr
 453         do ilam2=-180,179,iincr
 454           do ilam3=-180,179,iincr
 455             do ilam4=-180,179,iincr
 456               lambda1=ilam1*conv+delta2
 457               lambda2=ilam2*conv+delta2
 458               lambda3=ilam3*conv+delta2
 459               lambda4=ilam4*conv+delta2
 460               ele1=ele(lambda1,si1*lambda4,a1)
 461               eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
 462               eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
 463               if (si1.gt.0) then
 464                 eloc3=elocal(ity3,lambda3,gamma3-pi-lambda4,.false.)
 465               else
 466                 eloc3=elocal(ity3,lambda3,lambda4,.false.)
 467               endif
 468               sum1=sum1+ele1*eloc1*eloc2*eloc3
 469             enddo
 470           enddo
 471         enddo
 472       enddo
 473       pom=-1.0d0/(2*pi)**4*delta**4
 474       ene1=sum1*pom
 475       return
 476       end
 477 c-------------------------------------------------------------------------
 478       double precision function elocal(i,x,y,transp)
 479       implicit real*8 (a-h,o-z)
 480       include 'DIMENSIONS'
 481       include 'COMMON.TORSION'
 482       integer i
 483       double precision x,y,u(2),v(2),cu(2),dv(2),ev(2)
 484       double precision scalar2
 485       logical transp
 486       u(1)=dcos(x)
 487       u(2)=dsin(x)
 488       v(1)=dcos(y)
 489       v(2)=dsin(y)
 490       if (transp) then
 491         call matvec2(cc(1,1,i),v,cu)
 492         call matvec2(dd(1,1,i),u,dv)
 493         call matvec2(ee(1,1,i),u,ev)
 494         elocal=scalar2(b1(1,i),v)+scalar2(b2(1,i),u)+scalar2(cu,v)+
 495      &   scalar2(dv,u)+scalar2(ev,v)
 496       else
 497         call matvec2(cc(1,1,i),u,cu)
 498         call matvec2(dd(1,1,i),v,dv)
 499         call matvec2(ee(1,1,i),v,ev)
 500         elocal=scalar2(b1(1,i),u)+scalar2(b2(1,i),v)+scalar2(cu,u)+
 501      &   scalar2(dv,v)+scalar2(ev,u)
 502       endif
 503       return
 504       end
 505 c-------------------------------------------------------------------------
 506       double precision function ele(x,y,a)
 507       implicit none
 508       double precision x,y,a(2,2),si1,si2,u(2),v(2),av(2)
 509       double precision scalar2
 510       u(1)=-cos(x)
 511       u(2)= sin(x)
 512       v(1)=-cos(y)
 513       v(2)= sin(y)
 514       call matvec2(a,v,av)
 515       ele=scalar2(u,av)
 516       return
 517       end