pro change_cor,Ra0,dec0,Ra1,Dec1,
rapol,decpol,PA,j
        ; coordinate change betweeb (Ra0,Dec0) and (Ra1,Dec1)
        ; the pole of coord1 in coord0 is (pole_Ra0,pole_Dec0)
        ; the position angle of coord1 is PA
        ; generalized version  from the conversion from Galactic to
celestial
        ; j= 1 coord0 --> coord1
        ; j=2 coord1 --> coord0
        ; all in degrees


        radeg = 180.0d/!DPI
        sdp=sin(decpol/radeg)
        cdp=sqrt(1.0d0-sdp*sdp)

case j of
    1:  begin
                ras = ra0/radeg - rapol/radeg
                   sdec = sin(dec0/radeg)
                cdec = sqrt(1.0d0-sdec*sdec)
            sgb = sdec*sdp + cdec*cdp*cos(ras)

                Dec1 = radeg * asin(sgb)
            cgb = sqrt(1.0d0-sgb*sgb)
            sine = cdec * sin(ras) / cgb
            cose = (sdec-sdp*sgb) / (cdp*cgb)
            Ra1 = PA - radeg*atan(sine,cose)
            ltzero=where(Ra1 lt 0.0, Nltzero)
            if Nltzero ge 1 then Ra1[ltzero]=Ra1[ltzero]+360.0d0
                 return
    end
        2: begin
                sgb = sin(Dec1/radeg)
        cgb = sqrt(1.0d0-sgb*sgb)
        sdec = sgb*sdp + cgb*cdp*cos((PA-Ra1)/radeg)
        dec0 = radeg * asin(sdec)
        cdec = sqrt(1.0d0-sdec*sdec)
        sinf = cgb * sin((PA-RA1)/radeg) / cdec
        cosf = (sgb-sdp*sdec) / (cdp*cdec)
        ra0 = rapol + radeg*atan(sinf,cosf)
            return
        end
endcase
end